Prognostic and Diagnostic value of malaria pigment
James Watson
12/11/2021
1 Systematic review
sysreview=read.csv('~/Dropbox/Apps/Overleaf/Pigment_Prognosis_Severe_Malaria/SystematicReview/All_Identified_papers_merged.csv')
table(sysreview$Screening.pass)##
## No Yes
## 183 19table(sysreview$ReasonExclude)##
## Case report Chickens
## 19 12 1
## Commentary Culture Duplicated
## 1 33 3
## Ex vivo Fatal cases only In vitro
## 2 6 15
## Mice Monkey Not falciparum
## 35 3 7
## Not severe malaria Pigment not measured Puffins
## 25 12 1
## Rabbit Retrospective Review
## 1 2 23
## Subcutanaeous not blood
## 1sum(sysreview$N..Severe.malaria.[sysreview$Screening.pass=='Yes'])## [1] 27958median(sysreview$N..Severe.malaria.[sysreview$Screening.pass=='Yes'])## [1] 67range(sysreview$N..Severe.malaria.[sysreview$Screening.pass=='Yes'])## [1] 22 262962 Data
load('RData/Merged_pigment_data.RData')
pigmt_data$meta_code = apply(pigmt_data[, c('study','country')], 1,
function(x) paste(x[1],x[2],sep=': '))
table(pigmt_data$study)##
## AQ AQUAMAT Lyke et al SEAQUAMAT SMAC
## 560 5425 172 1461 26235# remove patients without an outcome variable
pigmt_data = filter(pigmt_data, !is.na(outcome))
my_cols_AQ = brewer.pal(n=9, name = 'Paired')
my_cols_SQ = brewer.pal(n=4, name = 'Accent')
my_cols_SMAC = brewer.pal(n=6, name = 'Set1')
aqViet_ind = pigmt_data$study == 'AQ'
aq_ind = pigmt_data$study == 'AQUAMAT'
sq_ind = pigmt_data$study == 'SEAQUAMAT'
smac_ind = pigmt_data$study == 'SMAC'
lyke_ind = pigmt_data$study == 'Lyke et al'
pigmt_data$site_f = as.factor(pigmt_data$site)
pigmt_data$country_f = as.factor(pigmt_data$country)
pigmt_data$study_f = as.factor(pigmt_data$study)
pigmt_data$prop_pigneut = pigmt_data$pig_neut/pigmt_data$pig_neut_denom
pigmt_data$prop_pigmono = pigmt_data$pig_mono/pigmt_data$pig_mono_denom
writeLines('Available PMN data:')## Available PMN data:table(pigmt_data$study, !is.na(pigmt_data$pig_neut))##
## FALSE TRUE
## AQ 77 483
## AQUAMAT 1214 4211
## Lyke et al 6 166
## SEAQUAMAT 128 1333
## SMAC 96 25845writeLines('Available PMM data:')## Available PMM data:table(pigmt_data$study, !is.na(pigmt_data$pig_mono))##
## FALSE TRUE
## AQ 259 301
## AQUAMAT 1239 4186
## Lyke et al 6 166
## SEAQUAMAT 129 1332
## SMAC 916 25025writeLines('Available parasite count data:')## Available parasite count data:table(pigmt_data$study, !is.na(pigmt_data$parasitaemia))##
## FALSE TRUE
## AQ 0 560
## AQUAMAT 639 4786
## Lyke et al 0 172
## SEAQUAMAT 0 1461
## SMAC 31 25910parameters for Box-Cox transformation
lambda_para = 50
lambda_pigment = 0.5/100
pigmt_data$lambda_para=lambda_para
pigmt_data$lambda_pigment=lambda_pigment2.1 Assess bias in data availibility
Availability of pigment data: is there a bias related to death? And related to parasitaemia?
mod_death_AQ=
glmer(outcome ~ as.numeric(!is.na(pig_neut))+(1|site:country_f),
family='binomial',
data = pigmt_data[aq_ind,])
summary(mod_death_AQ)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3042322 0.1674649 -13.759492 4.466553e-43
## as.numeric(!is.na(pig_neut)) -0.1292729 0.1145053 -1.128969 2.589110e-01mod_death_SQ=
glmer(outcome ~ as.numeric(!is.na(pig_neut))+(1|site:country_f),
family='binomial',
data = pigmt_data[sq_ind,])
summary(mod_death_SQ)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.9418980 0.3183736 -6.0994322 1.064459e-09
## as.numeric(!is.na(pig_neut)) 0.2560062 0.2653230 0.9648849 3.346025e-01mod_para_AQ =
glmer(as.numeric(!is.na(pig_neut)) ~ log10(parasitaemia+lambda_para)+
(1|site:country_f),
family='binomial',
data = pigmt_data[aq_ind,])
summary(mod_para_AQ)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.9576671 0.40784154 -4.800068 1.586119e-06
## log10(parasitaemia + lambda_para) 0.8245222 0.04590167 17.962794 3.811336e-72mod_para_SQ =
glm(as.numeric(!is.na(pig_neut)) ~ log10(parasitaemia+lambda_para) +
country_f,
family='binomial',
data = pigmt_data[sq_ind,])
summary(mod_para_SQ)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6297328 0.36391621 1.730433 8.355287e-02
## log10(parasitaemia + lambda_para) 0.3717608 0.07815339 4.756809 1.966770e-06
## country_fIndia 0.5375038 0.40031095 1.342716 1.793640e-01
## country_fIndonesia -0.7583936 0.22341400 -3.394566 6.873742e-04
## country_fMyanmar 1.0715903 0.26961268 3.974555 7.051107e-05mod_para =
glmer(as.numeric(!is.na(pig_neut)) ~ log10(parasitaemia+lambda_para) +
(1|site:country_f:study),
family='binomial',
data = pigmt_data[sq_ind|aq_ind,])
summary(mod_para)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5697068 0.36384541 -1.565794 1.173969e-01
## log10(parasitaemia + lambda_para) 0.7025905 0.03897889 18.024900 1.242372e-72par(las = 1, family='serif', mfrow=c(1,2))
hist(log10(pigmt_data$parasitaemia+lambda_para)[sq_ind], freq = F,
breaks = seq(0,7,by=.5), main = 'SEAQUAMAT', xaxt='n', yaxt='n',
xlab = 'Parasitaemia on admission', ylab='',ylim=c(0,.6))
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
ind_na = is.na(pigmt_data$pig_neut)
lines(density(log10(pigmt_data$parasitaemia+lambda_para)[ind_na&sq_ind],
na.rm = T),lwd=3,col='blue')
lines(density(log10(pigmt_data$parasitaemia+lambda_para)[!ind_na&sq_ind],
na.rm = T),lwd=3, col='red',lty=2)
legend('topleft', legend = c('Readable','Not readable'),
inset = 0.03, col = c('red','blue'),lwd=3,lty=c(2,1))
hist(log10(pigmt_data$parasitaemia+lambda_para)[aq_ind], freq = F,
breaks = seq(0,7,by=.5), main = 'AQUAMAT', xaxt='n', yaxt='n',
xlab = 'Parasitaemia on admission', ylab='',ylim=c(0,.6))
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
ind_na = is.na(pigmt_data$pig_neut)
lines(density(log10(pigmt_data$parasitaemia+lambda_para)[ind_na&aq_ind],
na.rm = T),lwd=3,col='blue')
lines(density(log10(pigmt_data$parasitaemia+lambda_para)[!ind_na&aq_ind],
na.rm = T),lwd=3, col='red',lty=2)
legend('topleft', legend = c('Readable','Not readable'),
inset = 0.03, col = c('red','blue'),lwd=3,lty=c(2,1))
table(high_para=pigmt_data$parasitaemia[aq_ind|sq_ind] > 10000,
available=!is.na(pigmt_data$pig_neut[aq_ind|sq_ind]))## available
## high_para FALSE TRUE
## FALSE 411 1010
## TRUE 540 42862.2 Overall summaries of data
Data summaries
xx=aggregate(outcome ~ meta_code, pigmt_data, mean)
xx$outcome = round(xx$outcome*100,1)
print(xx)## meta_code outcome
## 1 AQ: Vietnam 14.8
## 2 AQUAMAT: DRC 4.3
## 3 AQUAMAT: Ghana 4.8
## 4 AQUAMAT: Kenya 10.0
## 5 AQUAMAT: Mozambique 11.3
## 6 AQUAMAT: Nigeria 9.1
## 7 AQUAMAT: Rwanda 5.2
## 8 AQUAMAT: Tanzania 14.4
## 9 AQUAMAT: The Gambia 11.0
## 10 AQUAMAT: Uganda 6.3
## 11 Lyke et al: Mali 8.7
## 12 SEAQUAMAT: Bangladesh 28.0
## 13 SEAQUAMAT: India 23.9
## 14 SEAQUAMAT: Indonesia 9.3
## 15 SEAQUAMAT: Myanmar 14.4
## 16 SMAC: Gabon (Lambarene) 1.4
## 17 SMAC: Gabon (Libreville) 5.1
## 18 SMAC: Ghana 4.6
## 19 SMAC: Kenya 3.6
## 20 SMAC: Malawi 2.6
## 21 SMAC: The Gambia 9.4xx=aggregate(age ~ study, pigmt_data, quantile, probs=c(0,.5,1))
xx$age = round(xx$age,1)
rownames(xx$age)=xx$study
print(xx$age)## 0% 50% 100%
## AQ 15.0 30.0 79.0
## AQUAMAT 0.0 2.0 14.0
## Lyke et al 0.2 2.5 10.4
## SEAQUAMAT 2.0 25.0 87.0
## SMAC 0.0 2.2 15.0xx=aggregate(coma ~ study, pigmt_data, mean)
xx$coma = round(xx$coma*100)
print(xx)## study coma
## 1 AQ 52
## 2 AQUAMAT 36
## 3 Lyke et al 21
## 4 SEAQUAMAT 43
## 5 SMAC 9xx=aggregate(prop_pigneut ~ study, pigmt_data, function(x) sum(!is.na(x)))
print(xx)## study prop_pigneut
## 1 AQ 483
## 2 AQUAMAT 4211
## 3 Lyke et al 166
## 4 SEAQUAMAT 1333
## 5 SMAC 25845xx=aggregate(prop_pigmono ~ study, pigmt_data, function(x) sum(!is.na(x)))
print(xx)## study prop_pigmono
## 1 AQ 301
## 2 AQUAMAT 4186
## 3 Lyke et al 166
## 4 SEAQUAMAT 1332
## 5 SMAC 25025xx=aggregate(parasitaemia ~ study, pigmt_data, function(x) sum(!is.na(x)))
print(xx)## study parasitaemia
## 1 AQ 560
## 2 AQUAMAT 4786
## 3 Lyke et al 172
## 4 SEAQUAMAT 1461
## 5 SMAC 25910xx=aggregate(outcome ~ meta_code, pigmt_data[!is.na(pigmt_data$pig_neut),], mean)
xx$outcome = round(xx$outcome*100,1)
print(xx)## meta_code outcome
## 1 AQ: Vietnam 13.7
## 2 AQUAMAT: DRC 4.0
## 3 AQUAMAT: Ghana 4.1
## 4 AQUAMAT: Kenya 11.9
## 5 AQUAMAT: Mozambique 9.5
## 6 AQUAMAT: Nigeria 10.9
## 7 AQUAMAT: Rwanda 5.9
## 8 AQUAMAT: Tanzania 14.0
## 9 AQUAMAT: The Gambia 10.5
## 10 AQUAMAT: Uganda 5.6
## 11 Lyke et al: Mali 9.0
## 12 SEAQUAMAT: Bangladesh 27.8
## 13 SEAQUAMAT: India 23.9
## 14 SEAQUAMAT: Indonesia 10.2
## 15 SEAQUAMAT: Myanmar 14.8
## 16 SMAC: Gabon (Lambarene) 1.4
## 17 SMAC: Gabon (Libreville) 5.1
## 18 SMAC: Ghana 4.6
## 19 SMAC: Kenya 3.6
## 20 SMAC: Malawi 2.6
## 21 SMAC: The Gambia 9.5writeLines(sprintf('We have data on pigmented neutrophil counts in %s patients', sum(!is.na(pigmt_data$pig_neut))))## We have data on pigmented neutrophil counts in 32038 patientswriteLines(sprintf('We have data on pigmented monocyte counts in %s patients', sum(!is.na(pigmt_data$pig_mono))))## We have data on pigmented monocyte counts in 31010 patientswriteLines('Median pigmented neutrophil count by country and study:')## Median pigmented neutrophil count by country and study:aggregate(pig_neut ~ meta_code, data = pigmt_data, FUN = median)## meta_code pig_neut
## 1 AQ: Vietnam 2
## 2 AQUAMAT: DRC 2
## 3 AQUAMAT: Ghana 3
## 4 AQUAMAT: Kenya 4
## 5 AQUAMAT: Mozambique 2
## 6 AQUAMAT: Nigeria 2
## 7 AQUAMAT: Rwanda 2
## 8 AQUAMAT: Tanzania 3
## 9 AQUAMAT: The Gambia 2
## 10 AQUAMAT: Uganda 1
## 11 Lyke et al: Mali 3
## 12 SEAQUAMAT: Bangladesh 4
## 13 SEAQUAMAT: India 3
## 14 SEAQUAMAT: Indonesia 0
## 15 SEAQUAMAT: Myanmar 4
## 16 SMAC: Gabon (Lambarene) 1
## 17 SMAC: Gabon (Libreville) 2
## 18 SMAC: Ghana 0
## 19 SMAC: Kenya 0
## 20 SMAC: Malawi 0
## 21 SMAC: The Gambia 0writeLines('Median pigmented monocyte count by country and study:')## Median pigmented monocyte count by country and study:aggregate(pig_mono ~ meta_code, data = pigmt_data, FUN = median)## meta_code pig_mono
## 1 AQ: Vietnam 4.0
## 2 AQUAMAT: DRC 4.0
## 3 AQUAMAT: Ghana 3.0
## 4 AQUAMAT: Kenya 3.0
## 5 AQUAMAT: Mozambique 1.0
## 6 AQUAMAT: Nigeria 2.0
## 7 AQUAMAT: Rwanda 1.0
## 8 AQUAMAT: Tanzania 2.0
## 9 AQUAMAT: The Gambia 2.0
## 10 AQUAMAT: Uganda 1.0
## 11 Lyke et al: Mali 3.0
## 12 SEAQUAMAT: Bangladesh 1.0
## 13 SEAQUAMAT: India 1.0
## 14 SEAQUAMAT: Indonesia 0.0
## 15 SEAQUAMAT: Myanmar 1.0
## 16 SMAC: Gabon (Lambarene) 3.0
## 17 SMAC: Gabon (Libreville) 4.5
## 18 SMAC: Ghana 4.0
## 19 SMAC: Kenya 1.0
## 20 SMAC: Malawi 1.0
## 21 SMAC: The Gambia 5.0Some data plots
3 Prognostic value of parasite counts, PNs and PMs in AQUAMAT and SEAQUAMAT
3.1 Overall prognostic value
We use generalised additive models as the effects may be non-linear
## [1] 0.02994212 0.00000000## [1] 0 0## [1] 1.023784e-06 0.000000e+00## [1] 0.000000e+00 4.172026e-06## [1] 0.02535501 0.00000000## [1] 1.100682e-06 0.000000e+00## Warning in predict.gam(mod_PCT_SQ, newdata = data.frame(parasitaemia = 10^xs, :
## factor levels 0 not in original fit## Warning in predict.gam(mod_PCT_AQ, newdata = data.frame(parasitaemia = 10^xs, :
## factor levels 0 not in original fit## Warning in predict.gam(mod_PN_SQ, newdata = data.frame(prop_pigneut = xs, :
## factor levels 0 not in original fit## Warning in predict.gam(mod_PN_AQ, newdata = data.frame(prop_pigneut = xs, :
## factor levels 0 not in original fit## Warning in predict.gam(mod_PM_SQ, newdata = data.frame(prop_pigmono = xs, :
## factor levels 0 not in original fit## Warning in predict.gam(mod_PM_AQ, newdata = data.frame(prop_pigmono = xs, :
## factor levels 0 not in original fit
3.1.1 5% prognostic value
## Neutrophils
mod_5percentMORU_PMN =
glmer(outcome ~ (prop_pigneut>0.05) + (1|country_f:study),
family='binomial', data = pigmt_data[pigmt_data$study %in% c('AQ','AQUAMAT','SEAQUAMAT'), ])
xx=summary(mod_5percentMORU_PMN)
writeLines(sprintf('In the MORU data (n=%s), the odds ratio for death in those with >5%% pigment containing PMNs is %s (95%% CI %s to %s)',
length(xx$residuals),
round(exp(xx$coefficients[2,1]),2),
round(exp(xx$coefficients[2,1]-1.96*xx$coefficients[2,2]),2),
round(exp(xx$coefficients[2,1]+1.96*xx$coefficients[2,2]),2)))## In the MORU data (n=6027), the odds ratio for death in those with >5% pigment containing PMNs is 2.39 (95% CI 2.02 to 2.81)## Monocytes
mod_5percentMORU_PMM =
glmer(outcome ~ (prop_pigmono>0.05) + (1|country_f:study),
family='binomial', data = pigmt_data[pigmt_data$study %in% c('AQ','AQUAMAT','SEAQUAMAT'), ])
xx=summary(mod_5percentMORU_PMM)
writeLines(sprintf('In the MORU data (n=%s), the odds ratio for death in those with >5%% pigment containing monocytes is %s (95%% CI %s to %s)',
length(xx$residuals),
round(exp(xx$coefficients[2,1]),2),
round(exp(xx$coefficients[2,1]-1.96*xx$coefficients[2,2]),2),
round(exp(xx$coefficients[2,1]+1.96*xx$coefficients[2,2]),2)))## In the MORU data (n=5819), the odds ratio for death in those with >5% pigment containing monocytes is 1.31 (95% CI 1.1 to 1.55)3.1.2 In cerebral malaria
ind_cerebral = pigmt_data$coma==1
# Parasite counts
mod_PCT_AQ = gam(outcome ~ s(log10(parasitaemia+lambda_para),k=4) +
s(country_f, bs='re'),
family='binomial',
data = pigmt_data[aq_ind&ind_cerebral, ])
summary(mod_PCT_AQ)##
## Family: binomial
## Link function: logit
##
## Formula:
## outcome ~ s(log10(parasitaemia + lambda_para), k = 4) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5661 0.1822 -8.595 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(parasitaemia + lambda_para)) 1.449 1.752 0.828 0.486
## s(country_f) 6.402 8.000 67.241 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0433 Deviance explained = 4.88%
## UBRE = -0.089896 Scale est. = 1 n = 1680mod_PCT_SQ = gam(outcome ~ s(log10(parasitaemia+lambda_para),k=4) +
s(country_f, bs='re'),
family='binomial',
data = pigmt_data[sq_ind&ind_cerebral, ])
summary(mod_PCT_SQ)##
## Family: binomial
## Link function: logit
##
## Formula:
## outcome ~ s(log10(parasitaemia + lambda_para), k = 4) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9338 0.2839 -3.289 0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(parasitaemia + lambda_para)) 2.564 2.857 27.98 1.87e-05 ***
## s(country_f) 2.582 3.000 23.14 6.70e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.084 Deviance explained = 7.42%
## UBRE = 0.18405 Scale est. = 1 n = 629## Neutrophils
mod_PN_AQ = gam(outcome ~ s(log10(prop_pigneut+lambda_pigment),k=4) +
s(country_f, bs='re'),
family='binomial',
data = pigmt_data[aq_ind&ind_cerebral, ])
summary(mod_PN_AQ)##
## Family: binomial
## Link function: logit
##
## Formula:
## outcome ~ s(log10(prop_pigneut + lambda_pigment), k = 4) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5252 0.2328 -6.551 5.7e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(prop_pigneut + lambda_pigment)) 1.000 1.001 9.321 0.00227 **
## s(country_f) 6.849 8.000 61.900 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0538 Deviance explained = 6.23%
## UBRE = -0.080885 Scale est. = 1 n = 1461mod_PN_SQ = gam(outcome ~ s(log10(prop_pigneut+lambda_pigment),k=4) +
s(country_f, bs='re'),
family='binomial',
data = pigmt_data[sq_ind&ind_cerebral, ])
summary(mod_PN_SQ)##
## Family: binomial
## Link function: logit
##
## Formula:
## outcome ~ s(log10(prop_pigneut + lambda_pigment), k = 4) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8471 0.1779 -4.761 1.93e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(prop_pigneut + lambda_pigment)) 1.0 1 29.83 < 2e-16 ***
## s(country_f) 1.9 3 10.83 0.000914 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0782 Deviance explained = 7.07%
## UBRE = 0.18196 Scale est. = 1 n = 580## Monocytes
mod_PM_AQ = gam(outcome ~ s(log10(prop_pigmono+lambda_pigment),k=4) + s(country_f, bs='re'),
family='binomial',
data = pigmt_data[aq_ind&ind_cerebral, ])
summary(mod_PM_AQ)##
## Family: binomial
## Link function: logit
##
## Formula:
## outcome ~ s(log10(prop_pigmono + lambda_pigment), k = 4) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5122 0.2357 -6.417 1.39e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(prop_pigmono + lambda_pigment)) 1.506 1.805 0.841 0.652
## s(country_f) 6.880 8.000 68.372 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.05 Deviance explained = 5.85%
## UBRE = -0.076964 Scale est. = 1 n = 1452mod_PM_SQ = gam(outcome ~ s(log10(prop_pigmono+lambda_pigment),k=4) +
s(country_f, bs='re'),
family='binomial',
data = pigmt_data[sq_ind&ind_cerebral, ])
summary(mod_PM_SQ)##
## Family: binomial
## Link function: logit
##
## Formula:
## outcome ~ s(log10(prop_pigmono + lambda_pigment), k = 4) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8682 0.2458 -3.533 0.000411 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(prop_pigmono + lambda_pigment)) 1.483 1.765 8.929 0.0231 *
## s(country_f) 2.423 3.000 18.264 5.37e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0422 Deviance explained = 4.17%
## UBRE = 0.22185 Scale est. = 1 n = 580par(mfrow=c(3,2), las=1, bty='n', cex.lab=1.3,
cex.axis=1.3, family='serif')
####********* Parasite counts *********#######
xs = seq(log10(lambda_para),6,length.out = 100)
my_preds=predict(mod_PCT_AQ,
newdata=data.frame(parasitaemia=10^xs,
country_f = 0,
lambda_para=lambda_para),
exclude = "s(country_f)",se.fit = T)## Warning in predict.gam(mod_PCT_AQ, newdata = data.frame(parasitaemia = 10^xs, :
## factor levels 0 not in original fitplot(xs,100*inv.logit(my_preds$fit), panel.first = grid(),
ylab = 'Mortality (%)', xaxt='n', ylim = c(9,35),
xlab='Parasite count (per uL)',type='l',lwd=3)
axis(1, at = 1:6, labels = c(10,
expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
polygon(x = c(xs, rev(xs)),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
mtext(text = paste('AQUAMAT (n=',
sum(!is.na(mod_PCT_AQ$model[,2])),')',
sep=''),line = 1.5,
side = 3, adj = 0, cex=1.2)
my_preds=predict(mod_PCT_SQ,
newdata=data.frame(parasitaemia=10^xs,
country_f = 0,
lambda_para=lambda_para),
exclude = "s(country_f)",se.fit = T)## Warning in predict.gam(mod_PCT_SQ, newdata = data.frame(parasitaemia = 10^xs, :
## factor levels 0 not in original fitplot(xs, 100*inv.logit(my_preds$fit), panel.first = grid(),
ylab = 'Mortality (%)', xaxt='n',ylim = c(10,60),
xlab='Parasite count (per uL)',type='l',lwd=3)
axis(1, at = 1:6, labels = c(10,
expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
polygon(x = c(xs, rev(xs)),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
mtext(text = paste('SEAQUAMAT (n=',
sum(!is.na(mod_PCT_SQ$model[,2])),')',
sep=''),line = 1.5,
side = 3, adj = 0, cex=1.2)
####********* Neutrophils *********#######
xs = (0: round(100*quantile(pigmt_data$prop_pigneut[aq_ind],
probs = .99, na.rm = T)))/100
my_preds=predict(mod_PN_AQ,
newdata=data.frame(prop_pigneut=xs,
lambda_pigment=lambda_pigment,
country_f = 0),
exclude = "s(country_f)",se.fit = T)## Warning in predict.gam(mod_PN_AQ, newdata = data.frame(prop_pigneut = xs, :
## factor levels 0 not in original fitplot(log10(xs+lambda_pigment),
100*inv.logit(my_preds$fit), panel.first = grid(),
ylab = 'Mortality (%)', xaxt='n', ylim = c(9,35),
xlab='Pigment containing PMNs (%)',type='l',lwd=3)
axis(1, at = log10(c(0,1,2,5,10)/100+lambda_pigment),
labels = c(0,1,2,5,10))
polygon(x = log10(c(xs, rev(xs))+lambda_pigment),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
mtext(text = paste('AQUAMAT (n=',
sum(!is.na(mod_PN_AQ$model[,2])),')',
sep=''),line = 1.5,
side = 3, adj = 0, cex=1.2)
xs = (0: round(100*quantile(pigmt_data$prop_pigneut[sq_ind],
probs = .99, na.rm = T)))/100
my_preds=predict(mod_PN_SQ,
newdata=data.frame(prop_pigneut=xs,
lambda_pigment=lambda_pigment,
country_f = 0),
exclude = "s(country_f)",se.fit = T)## Warning in predict.gam(mod_PN_SQ, newdata = data.frame(prop_pigneut = xs, :
## factor levels 0 not in original fitplot(log10(xs+lambda_pigment),
100*inv.logit(my_preds$fit), panel.first = grid(),
ylab = 'Mortality (%)', xaxt='n',ylim = c(10,60),
xlab='Pigment containing PMNs (%)',
type='l',lwd=3)
axis(1, at = log10(c(0,1,2,5,10,20)/100 +lambda_pigment),
labels = c(0,1,2,5,10,20))
polygon(x = log10(c(xs, rev(xs))+lambda_pigment),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
mtext(text = paste('SEAQUAMAT (n=',
sum(!is.na(mod_PN_SQ$model[,2])),')',
sep=''),line = 1.5,
side = 3, adj = 0, cex=1.2)
####********* Monocytes *********#######
xs = (0: round(100*quantile(pigmt_data$prop_pigmono[aq_ind],
probs = .99, na.rm = T)))/100
my_preds=predict(mod_PM_AQ,
newdata=data.frame(prop_pigmono=xs,
lambda_pigment=lambda_pigment,
country_f = 0),
exclude = "s(country_f)",se.fit = T)## Warning in predict.gam(mod_PM_AQ, newdata = data.frame(prop_pigmono = xs, :
## factor levels 0 not in original fitplot(log10(xs+lambda_pigment),
100*inv.logit(my_preds$fit), panel.first = grid(),
ylab = 'Mortality (%)', xaxt='n', ylim = c(9,35),
xlab = 'Pigment containing monocytes (%)',
type='l',lwd=3)
axis(1, at = log10(c(0,1,5,10,20,40)/100+lambda_pigment),
labels = c(0,1,5,10,20,40))
polygon(x = log10(c(xs, rev(xs))+lambda_pigment),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
mtext(text = paste('AQUAMAT (n=',
sum(!is.na(mod_PM_AQ$model[,2])),')',
sep=''),line = 1.2,
side = 3, adj = 0, cex=1.2)
xs = (0: round(100*quantile(pigmt_data$prop_pigmono[sq_ind],
probs = .99, na.rm = T)))/100
my_preds=predict(mod_PM_SQ,
newdata=data.frame(prop_pigmono=xs,
lambda_pigment=lambda_pigment,
country_f = 0),
exclude = "s(country_f)",se.fit = T)## Warning in predict.gam(mod_PM_SQ, newdata = data.frame(prop_pigmono = xs, :
## factor levels 0 not in original fitplot(log10(xs+lambda_pigment),
100*inv.logit(my_preds$fit), panel.first = grid(),
ylab = 'Mortality (%)', xaxt='n',ylim = c(10,60),
xlab='Pigment containing monocytes (%)',type='l',lwd=3)
axis(1, at = log10(c(0,1,2,5,10,15)/100+lambda_pigment),
labels = c(0,1,2,5,10,15))
polygon(x = log10(c(xs, rev(xs))+lambda_pigment),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
mtext(text = paste('SEAQUAMAT (n=',
sum(!is.na(mod_PM_SQ$model[,2])),')',
sep=''),line = 1.5,
side = 3, adj = 0, cex=1.2)
3.2 Overall model
3.2.1 PMNs
ind_not_na_neut = !is.na(pigmt_data$prop_pigneut)
ind_MORU = aqViet_ind | aq_ind | sq_ind
writeLines(sprintf('We have PMN data on a total of %s patients from the 3 MORU trials', sum(ind_MORU & ind_not_na_neut)))## We have PMN data on a total of 6027 patients from the 3 MORU trialsMixed effects for site, country and study
# box-cox transformation
mod_all=glmer(outcome ~ log10(prop_pigneut+lambda_pigment) +
(1|site:country:study),
family='binomial', data = pigmt_data[ind_MORU,])
summary(mod_all)$coefficients## Estimate Std. Error z value
## (Intercept) -0.7189118 0.17582412 -4.088812
## log10(prop_pigneut + lambda_pigment) 0.9267863 0.08753163 10.588016
## Pr(>|z|)
## (Intercept) 4.335884e-05
## log10(prop_pigneut + lambda_pigment) 3.386948e-26writeLines('Odds-ratio for death for a tenfold increase in PMNs (box-cox transformation):')## Odds-ratio for death for a tenfold increase in PMNs (box-cox transformation):round(exp(summary(mod_all)$coefficients[2,1]+c(-1,0,1)*1.96*summary(mod_all)$coefficients[2,2]),2)## [1] 2.13 2.53 3.00# box-cox transformation - just using >0 data
writeLines('Odds-ratio for death for a tenfold increase in PMNs (box-cox transformation only >0 counts):')## Odds-ratio for death for a tenfold increase in PMNs (box-cox transformation only >0 counts):mod_all_pos=glmer(outcome ~ log10(prop_pigneut) +
(1|site:country:study),
family='binomial',
data = pigmt_data[ind_MORU&pigmt_data$prop_pigneut>0, ])
exp(summary(mod_all_pos)$coefficients)[2,1]## [1] 2.689449writeLines('Odds-ratio for death for PMNs > 5%:')## Odds-ratio for death for PMNs > 5%:mod_all=glmer(outcome ~ as.numeric(prop_pigneut > 0.05) +
(1|site:country:study),
family='binomial', data = pigmt_data[ind_MORU,])
summary(mod_all)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3897853 0.13501306 -17.70040 4.163329e-70
## as.numeric(prop_pigneut > 0.05) 0.8750535 0.08444453 10.36246 3.673714e-25round(exp(summary(mod_all)$coefficients[2,1]+c(-1,0,1)*1.96*summary(mod_all)$coefficients[2,2]),2)## [1] 2.03 2.40 2.833.2.2 monocytes
ind_not_na_mono = !is.na(pigmt_data$prop_pigmono)
writeLines(sprintf('We have PMN data on a total of %s patients from the 3 MORU trials', sum(ind_MORU & ind_not_na_mono)))## We have PMN data on a total of 5819 patients from the 3 MORU trialsMixed effects for site, country and study
# box-cox transformation
mod_all=glmer(outcome ~ log10(prop_pigmono+lambda_pigment) +
(1|site:country:study),
family='binomial', data = pigmt_data[ind_MORU,])
summary(mod_all)$coefficients## Estimate Std. Error z value
## (Intercept) -1.6179278 0.17487642 -9.251835
## log10(prop_pigmono + lambda_pigment) 0.3309347 0.06866728 4.819394
## Pr(>|z|)
## (Intercept) 2.206715e-20
## log10(prop_pigmono + lambda_pigment) 1.439947e-06writeLines('Odds-ratio for death for a tenfold increase in monocytes (box-cox transformation):')## Odds-ratio for death for a tenfold increase in monocytes (box-cox transformation):round(exp(summary(mod_all)$coefficients[2,1]+c(-1,0,1)*1.96*summary(mod_all)$coefficients[2,2]),2)## [1] 1.22 1.39 1.59# box-cox transformation - just using >0 data
writeLines('Odds-ratio for death for a tenfold increase in monocytes (box-cox transformation only >0 counts):')## Odds-ratio for death for a tenfold increase in monocytes (box-cox transformation only >0 counts):mod_all_pos=glmer(outcome ~ log10(prop_pigmono) +
(1|site:country:study),
family='binomial',
data = pigmt_data[ind_MORU&pigmt_data$prop_pigmono>0, ])
exp(summary(mod_all_pos)$coefficients)[2,1]## [1] 1.630057writeLines('Odds-ratio for death for monocytes > 5%:')## Odds-ratio for death for monocytes > 5%:mod_all=glmer(outcome ~ as.numeric(prop_pigmono > 0.05) +
(1|site:country:study),
family='binomial', data = pigmt_data[ind_MORU,])
summary(mod_all)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2094892 0.1531057 -14.431138 3.295954e-47
## as.numeric(prop_pigmono > 0.05) 0.2715869 0.0877934 3.093478 1.978254e-03round(exp(summary(mod_all)$coefficients[2,1]+c(-1,0,1)*1.96*summary(mod_all)$coefficients[2,2]),2)## [1] 1.10 1.31 1.564 Meta-analysis
4.1 Prognostic value in African children
Af_ind = aq_ind|smac_ind|lyke_ind
writeLines(sprintf('We have data on pigment containing PMNs in %s African children',sum(Af_ind & !is.na(pigmt_data$prop_pigneut))))## We have data on pigment containing PMNs in 30222 African childrenwriteLines(sprintf('We have data on pigment containing monocytes in %s African children',sum(Af_ind & !is.na(pigmt_data$prop_pigmono))))## We have data on pigment containing monocytes in 29377 African childrenwriteLines(sprintf('We have data on parasite densities in %s African children',sum(Af_ind & !is.na(pigmt_data$parasitaemia))))## We have data on parasite densities in 30868 African children### Models using Parasite counts
mod_PCT = gam(outcome ~ s(log10(parasitaemia+lambda_para),k=4) +
s(country_f, bs='re')+s(study_f, bs='re'),
family='binomial',
data = pigmt_data[Af_ind, ])
summary(mod_PCT)$s.pv## [1] 0.0000000000 0.0000000000 0.0001050986## Neutrophils
mod_PN = gam(outcome ~ s(log10(prop_pigneut+lambda_pigment),k=4) +
s(country_f, bs='re')+s(study_f, bs='re'),
family='binomial', data = pigmt_data[Af_ind, ])
summary(mod_PN)$s.pv## [1] 0.000000 0.000000 0.482333## Monocytes
mod_PM = gam(outcome ~ s(log10(prop_pigmono+lambda_pigment),k=4) +
s(country_f, bs='re')+s(study_f, bs='re'),
family='binomial', data = pigmt_data[Af_ind, ])
summary(mod_PM)$s.pv## [1] 0.000000000 0.000000000 0.009103607raw counts/comparisons
pigmt_data$death = factor(pigmt_data$outcome, levels = 1:0)
table(pigmt_data$outcome[Af_ind])##
## 0 1
## 29860 1678writeLines('For pPMNs > 5percent:')## For pPMNs > 5percent:xx_pPMN_5 = table(pPMN_geq_5=pigmt_data$prop_pigneut[Af_ind]>=0.05,
death=pigmt_data$death[Af_ind])
print(aggregate(outcome ~ prop_pigneut>=0.05, data = pigmt_data[Af_ind,], mean))## prop_pigneut >= 0.05 outcome
## 1 FALSE 0.04430995
## 2 TRUE 0.10474537oddsRatio(xx_pPMN_5)## [1] 2.523505writeLines('For coma:')## For coma:xx_coma = table(coma=pigmt_data$coma[Af_ind],
death=pigmt_data$death[Af_ind])
print(aggregate(outcome ~ coma, data = pigmt_data[Af_ind,], function(x) c(round(100*mean(x),1), length(x))))## coma outcome.1 outcome.2
## 1 0 3.0 27115.0
## 2 1 19.3 4410.0oddsRatio(xx_coma)## [1] 7.619706writeLines('For acidosis (deep breathing - Kussmauls/respiratory distress)')## For acidosis (deep breathing - Kussmauls/respiratory distress)xx_acidosis = table(acidosis=pigmt_data$acidosis[Af_ind],
death=pigmt_data$death[Af_ind])
print(aggregate(outcome ~ acidosis, data = pigmt_data[Af_ind,], mean))## acidosis outcome
## 1 0 0.03440035
## 2 1 0.19013333oddsRatio(xx_acidosis)## [1] 6.589903pigmt_data$lab_acidosis = ifelse(test = is.na(pigmt_data$lactate),
yes = pigmt_data$bun>10,
no = pigmt_data$lactate>5)
xx_labacidosis = table(acidosis=pigmt_data$lab_acidosis[Af_ind],
death=pigmt_data$death[Af_ind])
print(aggregate(outcome ~ lab_acidosis, data = pigmt_data[Af_ind,], mean))## lab_acidosis outcome
## 1 FALSE 0.02464469
## 2 TRUE 0.11064061oddsRatio(xx_labacidosis)## [1] 4.923531writeLines('For haemoglobin < 5g/dL')## For haemoglobin < 5g/dLxx_hb_5 = table(lactate=pigmt_data$hb[Af_ind]<5,
death=pigmt_data$death[Af_ind])
print(aggregate(outcome ~ hb<5, data = pigmt_data[Af_ind,], mean))## hb < 5 outcome
## 1 FALSE 0.04461987
## 2 TRUE 0.07656747oddsRatio(xx_hb_5)## [1] 1.775363writeLines('For parasitaemia > 10,000 per uL')## For parasitaemia > 10,000 per uLxx_para_100 = table(lactate=pigmt_data$parasitaemia[Af_ind]>100000,
death=pigmt_data$death[Af_ind])
print(aggregate(outcome ~ parasitaemia>100000, data = pigmt_data[Af_ind,], mean))## parasitaemia > 1e+05 outcome
## 1 FALSE 0.05165813
## 2 TRUE 0.05234115oddsRatio(xx_para_100)## [1] 1.013952writeLines('For hypoglycaemia')## For hypoglycaemiaxx_hypo = table(hypo=pigmt_data$hypoglycaemia[Af_ind],
death=pigmt_data$death[Af_ind])
print(aggregate(outcome ~ hypoglycaemia, data = pigmt_data[Af_ind,], mean))## hypoglycaemia outcome
## 1 0 0.04143789
## 2 1 0.19235563oddsRatio(xx_hypo)## [1] 5.509438fixef(glmer(outcome ~ coma + (1|site), family = 'binomial', data = pigmt_data[Af_ind,]))## (Intercept) coma
## -3.517913 2.087549fixef(glmer(outcome ~ acidosis + (1|site), family = 'binomial', data = pigmt_data[Af_ind,]))## (Intercept) acidosis
## -3.116600 1.845477fixef(glmer(outcome ~ coma+acidosis+hypoglycaemia+ log2(prop_pigneut+0.005) + (1|site), family = 'binomial', data = pigmt_data[Af_ind,]))## (Intercept) coma
## -3.29888949 1.69872342
## acidosis hypoglycaemia
## 1.28200172 1.00785084
## log2(prop_pigneut + 0.005)
## 0.08912746par(las=1, bty='n', cex.lab=1.5, cex.axis=1.5, family='serif')
layout(mat = matrix(c(1,2,2,2,3,4,4,4,5,6,6,6),nrow = 4,byrow = F))
####********* Parasite counts *********#######
xs = seq(log10(lambda_para),6.5,length.out = 100)
par(mar = c(1,6,5,4))
hist(log10(pigmt_data$parasitaemia[Af_ind]+50),
main ='', xlab='',xaxt='n',xlim = range(xs)+c(-.1,.1),
breaks = 30, ylab ='', yaxt='n')
axis(2,at = c(0,1000,2000,3000))
mtext(text = paste0('(n=',sum(notna_para&Af_ind),')'),
line = 1, side = 3, adj = 1, cex=1.2)
mtext(text = 'a', line = 1,side = 3, adj = 0, cex=1.2)
par(mar = c(5,6,0,4))
my_preds=predict(mod_PCT,
newdata=data.frame(parasitaemia=10^xs,
country_f = 0,
study_f = 0,
lambda_para=lambda_para),
exclude = c("s(country_f)","study_f"),
se.fit = T)## Warning in predict.gam(mod_PCT, newdata = data.frame(parasitaemia = 10^xs, :
## factor levels 0 not in original fit
## Warning in predict.gam(mod_PCT, newdata = data.frame(parasitaemia = 10^xs, :
## factor levels 0 not in original fitplot(xs, 100*inv.logit(my_preds$fit), panel.first = grid(),
ylab = 'Mortality (%)', xaxt='n',ylim = c(3,20),
xlab='Parasite count (per uL)',type='l',lwd=3)
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
polygon(x = c(xs, rev(xs)),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
####********* Neutrophils *********#######
xs = (0: round(100*quantile(pigmt_data$prop_pigneut,
probs = .99, na.rm = T)))/100
par(mar = c(1,6,5,4))
hist(log10(pigmt_data$prop_pigneut[Af_ind]+lambda_pigment),
main ='', xlab='', xaxt='n',ylab ='',
xlim = range(log10(xs+lambda_pigment))+c(-.1,.1), breaks = 20)
mtext(text = paste0('(n=',sum(notna_neut&Af_ind),')'),
line = 1, side = 3, adj = 1, cex=1.2)
mtext(text = 'b', line = 1,side = 3, adj = 0, cex=1.2)
par(mar = c(5,6,0,4))
my_preds=predict(mod_PN,
newdata=data.frame(prop_pigneut=xs,
lambda_pigment=lambda_pigment,
country_f = 0,
study_f=0),
exclude = c("s(country_f)","study_f)"),
se.fit = T)## Warning in predict.gam(mod_PN, newdata = data.frame(prop_pigneut = xs,
## lambda_pigment = lambda_pigment, : factor levels 0 not in original fit## Warning in predict.gam(mod_PN, newdata = data.frame(prop_pigneut = xs,
## lambda_pigment = lambda_pigment, : factor levels 0 not in original fitplot(log10(xs+lambda_pigment),
100*inv.logit(my_preds$fit), panel.first = grid(NA,NULL),
ylab = 'Mortality (%)', xaxt='n',ylim = c(3,20),
xlab='Pigment containing PMNs (%)',
type='l',lwd=3)
abline(v=log10(my_v_ticks/100+lambda_pigment), lty='dotted', col='lightgray')
axis(1, at = log10(c(0,1,5,20)/100 +lambda_pigment),
labels = c(0,1,5,20))
polygon(x = log10(c(xs, rev(xs))+lambda_pigment),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
####********* Monocytes *********#######
xs = (0: round(100*quantile(pigmt_data$prop_pigmono,
probs = .99, na.rm = T)))/100
par(mar = c(1,6,5,4))
hist(log10(pigmt_data$prop_pigmono[Af_ind]+lambda_pigment),
main ='', xlab='', xaxt='n',ylab ='',yaxt='n',
xlim = range(log10(xs+lambda_pigment))+c(-.1,.1), breaks = 20)
mtext(text = paste0('(n=',sum(notna_mono&Af_ind),')'),
line = 1, side = 3, adj = 1, cex=1.2)
mtext(text = 'c', line = 1,side = 3, adj = 0, cex=1.2)
axis(2, at=c(0,3000,6000,9000))
par(mar = c(5,6,0,4))
my_preds=predict(mod_PM,
newdata=data.frame(prop_pigmono=xs,
lambda_pigment=lambda_pigment,
country_f = 0,
study_f = 0),
exclude = c("s(country_f)","study_f"),
se.fit = T)## Warning in predict.gam(mod_PM, newdata = data.frame(prop_pigmono = xs,
## lambda_pigment = lambda_pigment, : factor levels 0 not in original fit## Warning in predict.gam(mod_PM, newdata = data.frame(prop_pigmono = xs,
## lambda_pigment = lambda_pigment, : factor levels 0 not in original fitplot(log10(xs+lambda_pigment),
100*inv.logit(my_preds$fit), panel.first = grid(NA,NULL),
ylab = 'Mortality (%)', xaxt='n',ylim = c(3,20),
xlab='Pigment containing monocytes (%)',type='l',lwd=3)
abline(v=log10(my_v_ticks/100+lambda_pigment), lty='dotted', col='lightgray')
axis(1, at = log10(c(0,1,5,20)/100+lambda_pigment),
labels = c(0,1,5,20))
polygon(x = log10(c(xs, rev(xs))+lambda_pigment),
y = 100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
border = NA, col = adjustcolor('grey',alpha.f = .3))
4.2 Prognostic value of pigmented neutrophils alone: meta-analysis
Remove patients with no pigment data
pigmt_data = filter(pigmt_data,
!is.na(prop_pigmono) | !is.na(prop_pigneut))
apply(pigmt_data, 2, function(x) round(100*mean(is.na(x)),1))## study site country id
## 0.0 0.0 0.0 81.2
## age art outcome resp_rate
## 0.2 0.0 0.0 0.9
## lactate acute_renal_failure bun hb
## 22.3 98.5 84.0 2.4
## hypoglycaemia acidosis base_excess coma
## 0.7 2.4 42.4 0.0
## parasitaemia pig_neut pig_neut_denom pig_mono
## 0.8 0.0 0.0 3.2
## pig_mono_denom meta_code site_f country_f
## 0.0 0.0 0.0 0.0
## study_f prop_pigneut prop_pigmono lambda_para
## 0.0 0.0 3.2 0.0
## lambda_pigment death lab_acidosis
## 0.0 0.0 7.0aqViet_ind = pigmt_data$study == 'AQ'
aq_ind = pigmt_data$study == 'AQUAMAT'
sq_ind = pigmt_data$study == 'SEAQUAMAT'
smac_ind = pigmt_data$study == 'SMAC'
lyke_ind = pigmt_data$study == 'Lyke et al'
Af_ind = aq_ind|smac_ind|lyke_indquantile(round(100*pigmt_data$prop_pigneut),
seq(0,1,by=.1), na.rm = T)## 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
## 0 0 0 0 0 0 1 2 3 6 100mycuts = c(-1,0,1,2,5,100)
pigmt_data$pig_neut100_bins =
cut(round(100*pigmt_data$prop_pigneut),
breaks = mycuts, right = T,
include.lowest = T)
table(pigmt_data$pig_neut100_bins)##
## [-1,0] (0,1] (1,2] (2,5] (5,100]
## 18847 2610 3761 3275 3545mod_all=glmer(outcome ~ pig_neut100_bins + (1|site:country:study),
family='binomial', data = pigmt_data)
xxs = summary(mod_all)$coefficients[-1, ]
quantile(round(100*pigmt_data$prop_pigmono), seq(0,1,by=.1), na.rm = T)## 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
## 0 0 0 0 0 2 3 4 7 13 100pigmt_data$pig_mono30_bins =
cut(round(100*pigmt_data$prop_pigmono),
breaks = mycuts, include.lowest = T)
table(pigmt_data$pig_mono30_bins)##
## [-1,0] (0,1] (1,2] (2,5] (5,100]
## 12645 1752 3707 4667 8239mod_all=glmer(outcome ~ pig_mono30_bins + (1|site:country:study),
family='binomial', data = pigmt_data)
xxs2 = summary(mod_all)$coefficients[-1, ]
par(mfrow=c(1,2), mar=c(5,6,4,6), las=1,
bty='n',family='serif', cex.lab=1.3, cex.axis=1.3)
plot(c(0,xxs[,1]), 0:nrow(xxs), pch=16, xlim = c(-0.25, 1.2),
ylab='', yaxt='n',cex=1.5, #panel.first = grid(),
xlab='Odds-ratio for death', xaxt='n')
abline(v=log(c(1,2,3)),h=0:5, lty='dotted', col='lightgray')
abline(v=0, lty=2, lwd=2)
mtext(text = 'a',side = 3,adj = 0,cex = 1.3, line = 1.5)
mtext(text = '(pigment containing PMNs)',side = 3,adj = 1,cex = 1.3, line = 1.5)
axis(1, at = log(c(1,2,3)), labels = 1:3)
axis(2, at = 0:nrow(xxs),
labels = c('0%\n(reference)','(0-1]%','(1-2]%','(2-5]%','>5%'),tick = F)
axis(4, at = 0:nrow(xxs),tick = F,
labels = paste('n',table(pigmt_data$pig_neut100_bins),
sep='='))
for(i in 1:nrow(xxs)){
upci = xxs[i,1]+1.96*xxs[i,2]
loci = xxs[i,1]-1.96*xxs[i,2]
lines(c(loci, upci), c(i,i), lwd=3)
}
rm(xxs)
plot(c(0,xxs2[,1]), 0:nrow(xxs2), pch=16, xlim = c(-0.25, 1.2),
ylab='', yaxt='n', cex=1.5, #panel.first = grid(),
xlab='Odds-ratio for death', xaxt='n')
abline(v=0, lty=2, lwd=2)
abline(v=log(c(1,2,3)),h=0:5, lty='dotted', col='lightgray')
mtext(text = 'b',side = 3,adj = 0,cex = 1.3, line = 1.5)
mtext(text = '(pigment containing monocytes)',side = 3,adj = 1,cex = 1.3, line = 1.5)
axis(1, at = log(c(1,2,3)), labels = 1:3)
axis(2, at = 0:nrow(xxs2),
labels = c('0%\n(reference)','(0-1]%','(1-2]%','(2-5]%','>5%'),tick = F)
axis(4, at = 0:nrow(xxs2),tick = F,
labels = paste('n',table(pigmt_data$pig_mono30_bins),
sep='='))
for(i in 1:nrow(xxs2)){
upci = xxs2[i,1]+1.96*xxs2[i,2]
loci = xxs2[i,1]-1.96*xxs2[i,2]
lines(c(loci, upci), c(i,i), lwd=3)
}
ind_not_na = !is.na(pigmt_data$prop_pigneut) & !is.na(pigmt_data$outcome)
ind_pos = pigmt_data$prop_pigneut > 0.05
event_e = c()
n_e = c()
event_c = c()
n_c = c()
lab = c()
for(ss in unique(pigmt_data$meta_code)){
ind = pigmt_data$meta_code==ss
my_study = ss
event_e = c(event_e,
sum(pigmt_data$outcome[ind_not_na&ind&ind_pos]))
n_e = c(n_e, sum(ind_not_na&ind&ind_pos))
event_c = c(event_c, sum(pigmt_data$outcome[ind_not_na&ind& !ind_pos]))
n_c = c(n_c, sum(ind_not_na&ind& !ind_pos))
lab = c(lab, ss)
}
mod=metabin(event.e = event_e, n.e = n_e,event.c = event_c,
n.c = n_c, studlab = lab,method = mymethod,
sm = 'OR',
comb.random = T)
par(family='serif')
forest(mod, lab.e = '>5% pPMN', lab.c = '<5% pPMN',
rightcols=myrightcols,
test.overall.random=F, print.pval.Q=F)
Just looking at SMAC - what were they thinking???
ind_not_na = !is.na(pigmt_data$prop_pigneut) & !is.na(pigmt_data$outcome) & pigmt_data$study=='SMAC'
ind_pos = pigmt_data$prop_pigneut > 0.05
event_e = c()
n_e = c()
event_c = c()
n_c = c()
lab = c()
for(ss in unique(pigmt_data$meta_code[pigmt_data$study=='SMAC'])){
ind = pigmt_data$meta_code==ss
my_study = ss
event_e = c(event_e,
sum(pigmt_data$outcome[ind_not_na&ind&ind_pos]))
n_e = c(n_e, sum(ind_not_na&ind&ind_pos))
event_c = c(event_c, sum(pigmt_data$outcome[ind_not_na&ind& !ind_pos]))
n_c = c(n_c, sum(ind_not_na&ind& !ind_pos))
lab = c(lab, ss)
}
mod=metabin(event.e = event_e, n.e = n_e,event.c = event_c,
n.c = n_c, studlab = lab,method = mymethod,
sm = 'OR',
comb.random = T)
par(family='serif')
forest(mod, lab.e = '>5% pPMN', lab.c = '<5% pPMN',
rightcols=myrightcols,
test.overall.random=F, print.pval.Q=F)
Difference between children and adults? (which is mostly a difference between high and low transmission settings).
children
ind_not_na = !is.na(pigmt_data$prop_pigneut) &
!is.na(pigmt_data$outcome) & !is.na(pigmt_data$age) &
pigmt_data$age <=15
ind_pos = pigmt_data$prop_pigneut > 0.05
event_e = c()
n_e = c()
event_c = c()
n_c = c()
lab = c()
for(ss in unique(pigmt_data$meta_code[!is.na(pigmt_data$age) & pigmt_data$age <= 15])){
ind = pigmt_data$meta_code==ss
my_study = ss
event_e = c(event_e,
sum(pigmt_data$outcome[ind_not_na&ind&ind_pos]))
n_e = c(n_e, sum(ind_not_na&ind&ind_pos))
event_c = c(event_c, sum(pigmt_data$outcome[ind_not_na&ind& !ind_pos]))
n_c = c(n_c, sum(ind_not_na&ind& !ind_pos))
lab = c(lab, ss)
}
mod=metabin(event.e = event_e, n.e = n_e,event.c = event_c,
n.c = n_c, studlab = lab,
sm = 'OR',method = mymethod,
comb.random = T)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00310898 (tol = 0.002, component 1)par(family='serif')
forest(mod, lab.e = '>5% pPMN', lab.c = '<5% pPMN',
rightcols=myrightcols,
test.overall.random=F, print.pval.Q=F)
Adults
ind_not_na = !is.na(pigmt_data$prop_pigneut) &
!is.na(pigmt_data$outcome) & !is.na(pigmt_data$age) &
pigmt_data$age >15
ind_pos = pigmt_data$prop_pigneut > 0.05
event_e = c()
n_e = c()
event_c = c()
n_c = c()
lab = c()
for(ss in unique(pigmt_data$meta_code[!is.na(pigmt_data$age) & pigmt_data$age > 15])){
ind = pigmt_data$meta_code==ss
my_study = ss
event_e = c(event_e,
sum(pigmt_data$outcome[ind_not_na&ind&ind_pos]))
n_e = c(n_e, sum(ind_not_na&ind&ind_pos))
event_c = c(event_c, sum(pigmt_data$outcome[ind_not_na&ind& !ind_pos]))
n_c = c(n_c, sum(ind_not_na&ind& !ind_pos))
lab = c(lab, ss)
}
mod=metabin(event.e = event_e, n.e = n_e,event.c = event_c,
n.c = n_c, studlab = lab,
sm = 'OR',method = mymethod,
comb.random = T)
par(family='serif')
forest(mod, lab.e = '>5% pPMN', lab.c = '<5% pPMN',
rightcols=myrightcols,
test.overall.random=F, print.pval.Q=F)
table(pigmt_data$study, pigmt_data$age>15)##
## FALSE TRUE
## AQ 8 475
## AQUAMAT 4195 0
## Lyke et al 166 0
## SEAQUAMAT 206 1127
## SMAC 25808 0pigmt_data$Africa = as.numeric(Af_ind)
mod1=glmer(outcome ~ as.numeric(pig_neut > 5) + (1|meta_code),
family = 'binomial', data = pigmt_data)
exp(fixef(mod1)[2])## as.numeric(pig_neut > 5)
## 2.245303summary(glmer(outcome ~ as.numeric(pig_neut >5) * (age<15) + (1|meta_code),
family = 'binomial', data = pigmt_data))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: outcome ~ as.numeric(pig_neut > 5) * (age < 15) + (1 | meta_code)
## Data: pigmt_data
##
## AIC BIC logLik deviance df.resid
## 13193.1 13235.0 -6591.6 13183.1 31972
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8490 -0.2889 -0.2004 -0.1616 8.6139
##
## Random effects:
## Groups Name Variance Std.Dev.
## meta_code (Intercept) 0.3164 0.5625
## Number of obs: 31977, groups: meta_code, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2139 0.2268 -9.763 < 2e-16 ***
## as.numeric(pig_neut > 5) 1.1852 0.1339 8.850 < 2e-16 ***
## age < 15TRUE -0.6548 0.2358 -2.777 0.00549 **
## as.numeric(pig_neut > 5):age < 15TRUE -0.4535 0.1468 -3.088 0.00201 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(_>5 a<15TR
## as.nmr(_>5) -0.284
## age<15TRUE -0.819 0.274
## a.(_>5):<15 0.265 -0.909 -0.292summary(glmer(outcome ~ as.numeric(pig_neut >5) * Africa + (1|site),
family = 'binomial', data = pigmt_data))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: outcome ~ as.numeric(pig_neut > 5) * Africa + (1 | site)
## Data: pigmt_data
##
## AIC BIC logLik deviance df.resid
## 13214.3 13256.1 -6602.1 13204.3 32033
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8310 -0.2887 -0.2002 -0.1612 8.6667
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.3325 0.5766
## Number of obs: 32038, groups: site, 29
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1990 0.2103 -10.456 < 2e-16 ***
## as.numeric(pig_neut > 5) 1.1580 0.1326 8.730 < 2e-16 ***
## Africa -0.6946 0.2571 -2.702 0.00689 **
## as.numeric(pig_neut > 5):Africa -0.4253 0.1463 -2.907 0.00365 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) as.(_>5) Africa
## as.nmr(_>5) -0.337
## Africa -0.824 0.277
## as.n(_>5):A 0.306 -0.906 -0.285summary(glmer(outcome ~ as.numeric(pig_neut >5) + (1|site),
family = 'binomial', data = pigmt_data[pigmt_data$age<15,]))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: outcome ~ as.numeric(pig_neut > 5) + (1 | site)
## Data: pigmt_data[pigmt_data$age < 15, ]
##
## AIC BIC logLik deviance df.resid
## 11726.4 11751.3 -5860.2 11720.4 30330
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5178 -0.2565 -0.2003 -0.1614 8.7006
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.3731 0.6108
## Number of obs: 30333, groups: site, 26
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.85562 0.14303 -19.96 <2e-16 ***
## as.numeric(pig_neut > 5) 0.73323 0.06149 11.93 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## as.nmr(_>5) -0.1444.3 Prognostic value of pigmented monocytes alone: meta-analysis
ind_not_na = !is.na(pigmt_data$prop_pigmono) & !is.na(pigmt_data$outcome)
ind_pos = pigmt_data$prop_pigmono > 0.05
event_e = c()
n_e = c()
event_c = c()
n_c = c()
lab = c()
for(ss in unique(pigmt_data$meta_code)){
ind = pigmt_data$meta_code==ss
my_study = pigmt_data$study[ind][1]
event_e = c(event_e, sum(pigmt_data$outcome[ind_not_na&ind&ind_pos]))
n_e = c(n_e, sum(ind_not_na&ind&ind_pos))
event_c = c(event_c, sum(pigmt_data$outcome[ind_not_na&ind& !ind_pos]))
n_c = c(n_c, sum(ind_not_na&ind& !ind_pos))
lab = c(lab, ss)
}
mod=metabin(event.e = event_e,
n.e = n_e,
event.c = event_c,
n.c = n_c, studlab = lab,
sm = 'OR',method = mymethod,
comb.random = T)
par(family='serif')
forest(mod, lab.e = '>5% pPMM', lab.c = '<5% pPMM',
rightcols=myrightcols,
test.overall.random=F, print.pval.Q=F)
Mixed effect model
mod_all=glmer(outcome ~ as.numeric(prop_pigmono) +
(1|site:country:study),
family='binomial', data = pigmt_data)
summary(mod_all)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.521557 0.1523346 -16.552750 1.529524e-61
## as.numeric(prop_pigmono) 2.030378 0.2501151 8.117774 4.748122e-16mod_all=glmer(outcome ~ log10(prop_pigmono+lambda_pigment) +
(1|site:country:study),
family='binomial', data = pigmt_data)
summary(mod_all)$coefficients## Estimate Std. Error z value
## (Intercept) -1.9207005 0.16216315 -11.844248
## log10(prop_pigmono + lambda_pigment) 0.2938478 0.04420657 6.647152
## Pr(>|z|)
## (Intercept) 2.304767e-32
## log10(prop_pigmono + lambda_pigment) 2.988183e-115 Prognostic value of pigmented monocytes in addition to severity data
5.1 Available clinical and lab data
# clinical variables
writeLines('Coma:')## Coma:table(pigmt_data$study,is.na(pigmt_data$coma))##
## FALSE TRUE
## AQ 483 0
## AQUAMAT 4217 0
## Lyke et al 166 0
## SEAQUAMAT 1333 0
## SMAC 25834 13writeLines('Acidosis:')## Acidosis:table(pigmt_data$study,is.na(pigmt_data$acidosis))##
## FALSE TRUE
## AQ 0 483
## AQUAMAT 4217 0
## Lyke et al 0 166
## SEAQUAMAT 1198 135
## SMAC 25847 0writeLines('Respiratory rate:')## Respiratory rate:table(pigmt_data$study,is.na(pigmt_data$resp_rate))##
## FALSE TRUE
## AQ 482 1
## AQUAMAT 4192 25
## Lyke et al 0 166
## SEAQUAMAT 1330 3
## SMAC 25754 93# lab variables
writeLines('Lactate:')## Lactate:table(pigmt_data$study,is.na(pigmt_data$lactate))##
## FALSE TRUE
## AQ 479 4
## AQUAMAT 0 4217
## Lyke et al 0 166
## SEAQUAMAT 0 1333
## SMAC 24414 1433writeLines('Glucose:')## Glucose:table(pigmt_data$study,is.na(pigmt_data$hypoglycaemia))##
## FALSE TRUE
## AQ 483 0
## AQUAMAT 4217 0
## Lyke et al 0 166
## SEAQUAMAT 1273 60
## SMAC 25847 0writeLines('Bun:')## Bun:table(pigmt_data$study,is.na(pigmt_data$bun))##
## FALSE TRUE
## AQ 228 255
## AQUAMAT 3650 567
## Lyke et al 0 166
## SEAQUAMAT 1252 81
## SMAC 0 25847pigmt_data$PMNs5 = as.numeric(pigmt_data$prop_pigneut > 0.05)
pigmt_data$PMMs5 = as.numeric(pigmt_data$prop_pigmono > 0.05)5.2 Fit models
5.2.1 AQ VIETNAM
Clinical: coma and respiratory rate Lab: lactate and glucose
#****** Neutrophils *********
# Model of prognostic value without severity included
mod_AQV_neut_alone = glm(outcome ~ PMNs5,
family='binomial',
data = pigmt_data[aqViet_ind, ])
summary(mod_AQV_neut_alone)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.381228 0.1819454 -13.087595 3.876887e-39
## PMNs5 1.783391 0.2829732 6.302332 2.931997e-10# Model of prognostic value with clinical markers of severity included
mod_AQV_neut_clinical = glm(outcome ~ PMNs5 +
coma + (resp_rate>35),
family='binomial',
data = pigmt_data[aqViet_ind, ])
summary(mod_AQV_neut_clinical)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7660733 0.2778625 -9.954828 2.402374e-23
## PMNs5 1.8360385 0.2941241 6.242394 4.309233e-10
## coma 0.3293968 0.2913084 1.130749 2.581605e-01
## resp_rate > 35TRUE 0.6932931 0.3145070 2.204381 2.749759e-02# Model of prognostic value with clinical & lab markers of severity included
mod_AQV_neut_lab = glm(outcome ~ PMNs5 +
coma + (resp_rate>35) + (lactate>5) + hypoglycaemia,
family='binomial',
data = pigmt_data[aqViet_ind, ])
summary(mod_AQV_neut_lab)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.0460083 0.3048899 -9.9905179 1.677033e-23
## PMNs5 1.3972170 0.3169699 4.4080426 1.043091e-05
## coma 0.2144539 0.3062394 0.7002818 4.837514e-01
## resp_rate > 35TRUE 0.1988386 0.3396067 0.5854969 5.582137e-01
## lactate > 5TRUE 1.4511100 0.3184874 4.5562561 5.207345e-06
## hypoglycaemia 0.4672507 0.4658315 1.0030466 3.158384e-01#****** Monocytes *********
# Model of prgnostic value without severity included
mod_AQV_mono_alone = glm(outcome ~ PMMs5,
family='binomial',
data = pigmt_data[aqViet_ind, ])
summary(mod_AQV_mono_alone)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.283414 0.5880758 -5.583318 2.359725e-08
## PMMs5 1.695835 0.6151301 2.756872 5.835712e-03# Model of prognostic value with clinical markers of severity included
mod_AQV_mono_clinical = glm(outcome ~ PMMs5 +
coma + (resp_rate>35),
family='binomial',
data = pigmt_data[aqViet_ind, ])
summary(mod_AQV_mono_clinical)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.36817038 0.6336792 -5.31526093 1.065046e-07
## PMMs5 1.63213338 0.6179726 2.64110963 8.263497e-03
## coma 0.02725479 0.3519870 0.07743125 9.382805e-01
## resp_rate > 35TRUE 0.55509619 0.4286262 1.29505882 1.952999e-01# Model of prognostic value with clinical & lab markers of severity included
mod_AQV_mono_lab = glm(outcome ~ PMMs5 +
coma + (resp_rate>35) + (lactate>5) + hypoglycaemia,
family='binomial',
data = pigmt_data[aqViet_ind, ])
summary(mod_AQV_mono_lab)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.85368982 0.6870465 -5.6090672 2.034200e-08
## PMMs5 1.33606879 0.6386269 2.0920960 3.642993e-02
## coma 0.06488358 0.3769613 0.1721227 8.633411e-01
## resp_rate > 35TRUE -0.29751024 0.4732309 -0.6286788 5.295594e-01
## lactate > 5TRUE 1.92759645 0.4006777 4.8108404 1.502970e-06
## hypoglycaemia 0.48746448 0.6140252 0.7938834 4.272633e-015.2.2 SEAQUAMAT
Clinical: coma and acidosis Lab: BUN and glucose
#****** Neutrophils *********
# Model of prgnostoc value without severity included
mod_SQ_neut_alone = glmer(outcome ~ PMNs5 +
(1|site:country),
family='binomial',
data = pigmt_data[sq_ind, ])
summary(mod_SQ_neut_alone)$coefficients[2,]## Estimate Std. Error z value Pr(>|z|)
## 1.013016e+00 1.484940e-01 6.821930e+00 8.982527e-12# Model of prognostic value with clinical markers of severity included
mod_SQ_neut_clinical = glmer(outcome ~ PMNs5 +
coma + acidosis +
(1|site:country),
family='binomial',
data = pigmt_data[sq_ind, ])
summary(mod_SQ_neut_clinical)$coefficients[2,]## Estimate Std. Error z value Pr(>|z|)
## 0.6160359021 0.1756288861 3.5076001205 0.0004521682# Model of prognostic value with clinical & lab markers of severity included
mod_SQ_neut_lab = glmer(outcome ~ PMNs5 +
coma + acidosis + (bun>10) +
hypoglycaemia + (1|site:country),
family='binomial',
data = pigmt_data[sq_ind, ])
summary(mod_SQ_neut_lab)$coefficients[2,]## Estimate Std. Error z value Pr(>|z|)
## 0.585524466 0.179631179 3.259592617 0.001115723#****** Monocytes *********
# Model of prgnostoc value without severity included
mod_SQ_mono_alone = glmer(outcome ~ PMMs5 +
(1|site:country),
family='binomial',
data = pigmt_data[sq_ind, ])
summary(mod_SQ_mono_alone)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8475029 0.2018076 -9.154772 5.447269e-20
## PMMs5 0.5142849 0.1508698 3.408799 6.524947e-04# Model of prognostic value with clinical markers of severity included
mod_SQ_mono_clinical = glmer(outcome ~ PMMs5 +
coma + acidosis +
(1|site:country),
family='binomial',
data = pigmt_data[sq_ind, ])
summary(mod_SQ_mono_clinical)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.9898304 0.2710736 -14.718623 4.895236e-49
## PMMs5 0.2687599 0.1824519 1.473045 1.407388e-01
## coma 1.3340290 0.1862776 7.161511 7.979272e-13
## acidosis 2.1896671 0.2250487 9.729748 2.251501e-22# Model of prognostic value with clinical & lab markers of severity included
mod_SQ_mono_lab = glmer(outcome ~ PMMs5 +
coma + acidosis + (bun>10) +
hypoglycaemia + (1|site:country),
family='binomial',
data = pigmt_data[sq_ind, ])
summary(mod_SQ_mono_lab)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.9308469 0.5156461 -9.562464 1.149861e-21
## PMMs5 0.1963372 0.1859240 1.056008 2.909647e-01
## coma 1.2672222 0.1894807 6.687868 2.264448e-11
## acidosis 2.0610195 0.2281334 9.034275 1.650925e-19
## bun > 10TRUE 1.1274449 0.4783478 2.356956 1.842541e-02
## hypoglycaemia 1.2684312 0.5305895 2.390607 1.682053e-025.2.3 AQUAMAT
Clinical: coma and acidosis Lab: BUN and glucose
#****** Neutrophils *********
# Model of prognostic value without severity included
mod_AQ_neut_alone = glmer(outcome ~ PMNs5 + (1|site:country),
family='binomial', data = pigmt_data[aq_ind, ])
summary(mod_AQ_neut_alone)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.589313 0.1517512 -17.062882 2.803771e-65
## PMNs5 0.651347 0.1113824 5.847844 4.979851e-09# Model of prognostic value with clinical markers of severity included
mod_AQ_neut_clinical = glmer(outcome ~ PMNs5 +
coma + acidosis +
(1|site:country),
family='binomial',
data = pigmt_data[aq_ind, ])
summary(mod_AQ_neut_clinical)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.6626794 0.1970765 -18.585064 4.244936e-77
## PMNs5 0.4903675 0.1187539 4.129276 3.639077e-05
## coma 1.7412533 0.1289991 13.498180 1.602880e-41
## acidosis 1.1712318 0.1315174 8.905525 5.313383e-19# Model of prognostic value with clinical & lab markers of severity included
mod_AQ_neut_lab = glmer(outcome ~ PMNs5 +
coma + acidosis + (bun>10) +
hypoglycaemia + (1|site:country),
family='binomial',
data = pigmt_data[aq_ind, ])
summary(mod_AQ_neut_lab)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.3304912 0.2297951 -18.845008 3.229141e-79
## PMNs5 0.2668282 0.1328897 2.007892 4.465472e-02
## coma 1.6483588 0.1408761 11.700768 1.263047e-31
## acidosis 0.7797244 0.1486133 5.246665 1.548769e-07
## bun > 10TRUE 0.8619773 0.1640212 5.255280 1.477998e-07
## hypoglycaemia 1.2303752 0.1484527 8.287993 1.151751e-16#****** Monocytes *********
# Model of prognostic value without severity included
mod_AQ_mono_alone = glmer(outcome ~ PMMs5 + (1|site:country),
family='binomial', data = pigmt_data[aq_ind, ])
summary(mod_AQ_mono_alone)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.45405542 0.1702992 -14.4102598 4.460292e-47
## PMMs5 0.09368014 0.1097358 0.8536876 3.932781e-01# Model of prognostic value with clinical markers of severity included
mod_AQ_mono_clinical = glmer(outcome ~ PMMs5 +
coma + acidosis +
(1|site:country),
family='binomial',
data = pigmt_data[aq_ind, ])
summary(mod_AQ_mono_clinical)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.528956036 0.2058681 -17.14182654 7.234202e-66
## PMMs5 -0.009519482 0.1161335 -0.08197013 9.346705e-01
## coma 1.777158114 0.1289767 13.77890796 3.413974e-43
## acidosis 1.172779563 0.1314186 8.92400088 4.497403e-19# Model of prognostic value with clinical & lab markers of severity included
mod_AQ_mono_lab = glmer(outcome ~ PMMs5 +
coma + acidosis + (bun>10) +
hypoglycaemia + (1|site:country),
family='binomial',
data = pigmt_data[aq_ind, ])
summary(mod_AQ_mono_lab)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.22837171 0.2401853 -17.6046237 2.270111e-69
## PMMs5 -0.08535355 0.1298821 -0.6571617 5.110770e-01
## coma 1.66810435 0.1412783 11.8072195 3.582122e-32
## acidosis 0.76565470 0.1493604 5.1262244 2.956103e-07
## bun > 10TRUE 0.88365483 0.1640855 5.3853301 7.231180e-08
## hypoglycaemia 1.25272875 0.1486131 8.4294635 3.472537e-175.2.4 SMAC
Clinical: coma and acidosis (deep breathing) Lab: lactate and glucose
#**** Neutrophils *****
# Model of prognostic value without severity included
mod_SMAC_neut_alone = glmer(outcome ~ PMNs5 + (1|site:country),
family='binomial',
data = pigmt_data[smac_ind, ])
summary(mod_SMAC_neut_alone)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.3218817 0.25889509 -12.830995 1.099426e-37
## PMNs5 0.9249743 0.09287892 9.958926 2.305405e-23# Model of prognostic value with clinical markers of severity included
mod_SMAC_neut_clinical = glmer(outcome ~ PMNs5 +
coma + acidosis + (1|site:country),
family='binomial',
data = pigmt_data[smac_ind, ])
summary(mod_SMAC_neut_clinical)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.9417665 0.23445370 -16.81256 1.974612e-63
## PMNs5 0.4976095 0.10158591 4.89841 9.661506e-07
## coma 1.7411385 0.07154428 24.33652 8.054366e-131
## acidosis 1.4891664 0.07050895 21.12025 5.182394e-99# Model of prognostic value with clinical & lab markers of severity included
mod_SMAC_neut_lab = glmer(outcome ~ PMNs5 +
coma + acidosis + (lactate>5) +
hypoglycaemia + (1|site:country),
family='binomial',
data = pigmt_data[smac_ind, ])
summary(mod_SMAC_neut_lab)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.3303018 0.22319318 -19.401587 7.482620e-84
## PMNs5 0.2395052 0.10732836 2.231518 2.564683e-02
## coma 1.5508821 0.07540598 20.567098 5.410877e-94
## acidosis 1.1637290 0.07593542 15.325246 5.185766e-53
## lactate > 5TRUE 0.9417572 0.07335555 12.838255 1.001061e-37
## hypoglycaemia 0.8390839 0.08717468 9.625317 6.251403e-22#**** Monocytes *****
# Model of prognostic value without severity included
mod_SMAC_mono_alone = glmer(outcome ~ PMMs5 + (1|site:country),
family='binomial',
data = pigmt_data[smac_ind, ])
summary(mod_SMAC_mono_alone)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.2828836 0.23823730 -13.779889 3.367882e-43
## PMMs5 0.3134577 0.07030566 4.458498 8.253593e-06length(summary(mod_SMAC_mono_alone)$residuals)## [1] 25025# Model of prognostic value with clinical markers of severity included
mod_SMAC_mono_clinical = glmer(outcome ~ PMMs5 +
coma + acidosis +
(1|site:country),
family='binomial',
data = pigmt_data[smac_ind, ])
summary(mod_SMAC_mono_clinical)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.9190705 0.22382337 -17.509657 1.209201e-68
## PMMs5 0.1008064 0.07492402 1.345448 1.784807e-01
## coma 1.7542667 0.07271620 24.124838 1.371954e-128
## acidosis 1.4870576 0.07139577 20.828370 2.394748e-96length(summary(mod_SMAC_mono_clinical)$residuals)## [1] 25016# Model of prognostic value with clinical & lab markers of severity included
mod_SMAC_mono_lab = glmer(outcome ~ PMMs5 +
coma + acidosis + (lactate>5) +
hypoglycaemia + (1|site:country),
family='binomial',
data = pigmt_data[smac_ind, ])
summary(mod_SMAC_mono_lab)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.29372367 0.22259512 -19.289388 6.594870e-83
## PMMs5 -0.08685903 0.07998441 -1.085950 2.775013e-01
## coma 1.56417000 0.07677523 20.373368 2.881623e-92
## acidosis 1.14597953 0.07699265 14.884272 4.170033e-50
## lactate > 5TRUE 0.94631349 0.07463487 12.679240 7.707542e-37
## hypoglycaemia 0.84785613 0.08788049 9.647831 5.020645e-22length(summary(mod_SMAC_mono_lab)$residuals)## [1] 236515.3 Output
Make summaries
#**Neutrophils**
my_var = 'PMNs5'
coefs_non_adj = data.frame(
rbind(summary(mod_AQV_neut_alone)$coefficients[my_var, 1:2],
summary(mod_SQ_neut_alone)$coefficients[my_var, 1:2],
summary(mod_AQ_neut_alone)$coefficients[my_var, 1:2],
summary(mod_SMAC_neut_alone)$coefficients[my_var, 1:2]))
coefs_non_adj$upperCI = coefs_non_adj$Estimate+1.96*coefs_non_adj$Std..Error
coefs_non_adj$lowerCI = coefs_non_adj$Estimate-1.96*coefs_non_adj$Std..Error
coefs_adj = data.frame(
rbind(summary(mod_AQV_neut_clinical)$coefficients[my_var, 1:2],
summary(mod_SQ_neut_clinical)$coefficients[my_var, 1:2],
summary(mod_AQ_neut_clinical)$coefficients[my_var, 1:2],
summary(mod_SMAC_neut_clinical)$coefficients[my_var, 1:2]))
coefs_adj$upperCI = coefs_adj$Estimate+1.96*coefs_adj$Std..Error
coefs_adj$lowerCI = coefs_adj$Estimate-1.96*coefs_adj$Std..Error
coefs_adj2 = data.frame(
rbind(summary(mod_AQV_neut_lab)$coefficients[my_var, 1:2],
summary(mod_SQ_neut_lab)$coefficients[my_var, 1:2],
summary(mod_AQ_neut_lab)$coefficients[my_var, 1:2],
summary(mod_SMAC_neut_lab)$coefficients[my_var, 1:2]))
coefs_adj2$upperCI = coefs_adj2$Estimate+1.96*coefs_adj2$Std..Error
coefs_adj2$lowerCI = coefs_adj2$Estimate-1.96*coefs_adj2$Std..Error
#**Monocytes**
my_var = 'PMMs5'
coefs_non_adj_PM = data.frame(
rbind(summary(mod_AQV_mono_alone)$coefficients[my_var, 1:2],
summary(mod_SQ_mono_alone)$coefficients[my_var, 1:2],
summary(mod_AQ_mono_alone)$coefficients[my_var, 1:2],
summary(mod_SMAC_mono_alone)$coefficients[my_var, 1:2]))
coefs_non_adj_PM$upperCI = coefs_non_adj_PM$Estimate+1.96*coefs_non_adj_PM$Std..Error
coefs_non_adj_PM$lowerCI = coefs_non_adj_PM$Estimate-1.96*coefs_non_adj_PM$Std..Error
coefs_adj_PM = data.frame(
rbind(summary(mod_AQV_mono_clinical)$coefficients[my_var, 1:2],
summary(mod_SQ_mono_clinical)$coefficients[my_var, 1:2],
summary(mod_AQ_mono_clinical)$coefficients[my_var, 1:2],
summary(mod_SMAC_mono_clinical)$coefficients[my_var, 1:2]))
coefs_adj_PM$upperCI = coefs_adj_PM$Estimate+1.96*coefs_adj_PM$Std..Error
coefs_adj_PM$lowerCI = coefs_adj_PM$Estimate-1.96*coefs_adj_PM$Std..Error
coefs_adj2_PM = data.frame(
rbind(summary(mod_AQV_mono_lab)$coefficients[my_var, 1:2],
summary(mod_SQ_mono_lab)$coefficients[my_var, 1:2],
summary(mod_AQ_mono_lab)$coefficients[my_var, 1:2],
summary(mod_SMAC_mono_lab)$coefficients[my_var, 1:2]))
coefs_adj2_PM$upperCI = coefs_adj2_PM$Estimate+1.96*coefs_adj2_PM$Std..Error
coefs_adj2_PM$lowerCI = coefs_adj2_PM$Estimate-1.96*coefs_adj2_PM$Std..Errorseverity at baseline across studies
ind_not_na_neut = !is.na(pigmt_data$prop_pigneut)
ind_not_na_mono = !is.na(pigmt_data$prop_pigmono)
xlims = log(c(.7, 16))
my_cols = brewer.pal(n = 12, name = 'Paired')[c(2,10,6,8)]
par(mar=c(5,10,4,2), las=1, bty='n', mfrow=c(1,2),cex.lab=1.2, cex.axis=1.2)
epsilon = -0.1
plot(coefs_non_adj[,'Estimate'], (4:1)-epsilon, pch=17,
ylim = c(1+epsilon, 4-epsilon),cex=1.5,
xlim = xlims, ylab = '', yaxt='n', xaxt='n',col=my_cols[1],
xlab = expression('Odds-ratio for death'))
points(coefs_adj[,'Estimate'], (4:1), pch=16,cex=1.5,col=my_cols[2])
points(coefs_adj2[,'Estimate'], (4:1)+epsilon, pch=15,cex=1.5,col=my_cols[3])
axis(1, at = log(c(0.7,1,2,4,8,16)), labels = c(0.7,1,2,4,8,16))
axis(2, at = 4:1, labels = c(paste0('AQ Vietnam\nn=',sum(aqViet_ind&ind_not_na_neut)),
paste0('SEAQUAMAT\nn=',sum(sq_ind&ind_not_na_neut)),
paste0('AQUAMAT\nn=',sum(aq_ind&ind_not_na_neut)),
paste0('SMAC\nn=',sum(smac_ind&ind_not_na_neut))),tick = F)
abline(v=log(c(2,4,8,16)),h=0:5, lty='dotted', col='lightgray')
abline(v=0, lty=2, lwd=3)
for(j in 1:4){
i = rev(1:4)[j]
lines(c(coefs_non_adj$lowerCI[j], coefs_non_adj$upperCI[j]),
c(i-epsilon, i-epsilon), lwd=2,col=my_cols[1])
lines(c(coefs_adj$lowerCI[j], coefs_adj$upperCI[j]),
c(i, i), lwd=2,col=my_cols[2])
lines(c(coefs_adj2$lowerCI[j], coefs_adj2$upperCI[j]),
c(i+epsilon, i+epsilon), lwd=2,col=my_cols[3])
}
legend('right',
legend = c('Non-adjusted',
'Adjusted: clinical',
'Adjusted: clinical and lab'),
lwd = 3, col = my_cols, pch = 17:15, cex=1.2)
mtext(text = 'a (pigment containing PMNs >5%)',side = 3,adj = 0, cex = 1.3, line=1.5)
par(mar=c(5,10,4,2), las=1)
plot(coefs_non_adj_PM[,'Estimate'], (4:1)-epsilon, pch=17,
ylim = c(1+epsilon, 4-epsilon),cex=1.5,
xlim = xlims, ylab = '', yaxt='n', xaxt='n',col=my_cols[1],
xlab = expression('Odds-ratio for death'))
points(coefs_adj_PM[,'Estimate'], (4:1), cex=1.5,pch=16,col=my_cols[2])
points(coefs_adj2_PM[,'Estimate'], (4:1)+epsilon,cex=1.5, pch=15,col=my_cols[3])
abline(v=log(c(2,4,8,16)),h=0:5, lty='dotted', col='lightgray')
axis(1, at = log(c(0.7,1,2,4,8,16)), labels = c(0.7,1,2,4,8,16))
axis(2, at = 4:1, labels = c(paste0('AQ Vietnam\nn=',sum(aqViet_ind&ind_not_na_mono)),
paste0('SEAQUAMAT\nn=',sum(sq_ind&ind_not_na_mono)),
paste0('AQUAMAT\nn=',sum(aq_ind&ind_not_na_mono)),
paste0('SMAC\nn=',sum(smac_ind&ind_not_na_mono))),tick = F)
abline(v=0, lty=2, lwd=3)
for(j in 1:4){
i = rev(1:4)[j]
lines(c(coefs_non_adj_PM$lowerCI[j], coefs_non_adj_PM$upperCI[j]),
c(i-epsilon, i-epsilon), lwd=2,col=my_cols[1])
lines(c(coefs_adj_PM$lowerCI[j], coefs_adj_PM$upperCI[j]),
c(i, i), lwd=2,col=my_cols[2])
lines(c(coefs_adj2_PM$lowerCI[j], coefs_adj2_PM$upperCI[j]),
c(i+epsilon, i+epsilon), lwd=2,col=my_cols[3])
}
mtext(text = 'b (pigment containing monocytes >5%)',side = 3,adj = 0, cex = 1.3, line=1.5)
6 Parasitaemia and pigment
par(las=1, mfrow=c(2,2),family='serif',cex.lab=1.3, cex.axis=1.3)
plot(log10(pigmt_data$parasitaemia+lambda_para),
jitter(log10(pigmt_data$prop_pigneut+lambda_pigment),
amount = 0.02),ylab = 'Pigment containing PMNs (%)',
xlab = 'Parasite density (per uL)', xlim = c(1.5,6.2),
pch='.',xaxt='n', yaxt='n',panel.first=grid())
axis(2, at = log10(c(0,1,2,5,10,20,40)/100 +lambda_pigment),
labels = c(0,1,2,5,10,20,40))
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
abline(h = log10(0.05+lambda_pigment),lty=2,lwd=2, col='red')
mod = glmer(PMNs5 ~ log10(parasitaemia+lambda_para) + (1 |site:country:study), family='binomial', data = pigmt_data)
summary(mod)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.5289872 0.25391632 -21.77484 4.018742e-105
## log10(parasitaemia + lambda_para) 0.8853594 0.02725899 32.47953 2.074645e-231xs = 10^(seq(1.5, 6, length.out = 100))
preds=predict(mod, newdata=data.frame(parasitaemia=xs,
site = -1, country= -1, study= -1,
lambda_para=lambda_para),
type='response', allow.new.levels = T)
plot(log10(xs), preds, type='l',
ylab = 'Pigment containing PMNs >5% (probability)',
xlab = 'Parasite density (per uL)', xlim = c(1.5,6.2),
xaxt='n', panel.first=grid(),lwd=3)
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
#### Monocytes
plot(log10(pigmt_data$parasitaemia+lambda_para),
jitter(log10(pigmt_data$prop_pigmono+lambda_pigment),
amount = 0.02),ylab = 'Pigment containing monocytes (%)',
xlab = 'Parasite density (per uL)', xlim = c(1.5,6.2),
pch='.',xaxt='n', yaxt='n',panel.first=grid())
axis(2, at = log10(c(0,1,2,5,10,20,40)/100 +lambda_pigment),
labels = c(0,1,2,5,10,20,40))
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
abline(h = log10(0.05+lambda_pigment),lty=2,lwd=2, col='red')
mod = glmer(PMMs5 ~ log10(parasitaemia+lambda_para) + (1 |site:country:study), family='binomial', data = pigmt_data)
summary(mod)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7481587 0.18483609 -14.86808 5.311094e-50
## log10(parasitaemia + lambda_para) 0.5150772 0.01613894 31.91519 1.643639e-223xs = 10^(seq(1.5, 6, length.out = 100))
preds=predict(mod, newdata=data.frame(parasitaemia=xs,site = -1, country= -1, study= -1, lambda_para=lambda_para),type='response', allow.new.levels = T)
plot(log10(xs), preds, type='l',
ylab = 'Pigment containing monocytes >5% (probability)',
xlab = 'Parasite density (per uL)', xlim = c(1.5,6.2),
xaxt='n', panel.first=grid(),lwd=3)
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
par(mfrow=c(1,2))
mod = glmer(as.numeric(prop_pigneut>0) ~ log10(parasitaemia+lambda_para) + (1 |site:country:study), family='binomial', data = pigmt_data)
summary(mod)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3718346 0.24154917 -9.819262 9.302231e-23
## log10(parasitaemia + lambda_para) 0.7830075 0.01515935 51.651772 0.000000e+00preds=predict(mod, newdata=data.frame(parasitaemia=xs,site = -1, country= -1, study= -1, lambda_para=lambda_para),type='response', allow.new.levels = T)
plot(log10(xs), preds, type='l',
ylab = 'Pigment containing PMNs >0% (probability)',
xlab = 'Parasite density (per uL)', xlim = c(1.5,6.2),
xaxt='n', panel.first=grid(),lwd=3)
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
mod = glmer(as.numeric(prop_pigmono>0) ~ log10(parasitaemia+lambda_para) + (1 |site:country:study), family='binomial', data = pigmt_data)
summary(mod)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0457954 0.15190085 -13.46796 2.41429e-41
## log10(parasitaemia + lambda_para) 0.6426988 0.01349885 47.61138 0.00000e+00preds=predict(mod, newdata=data.frame(parasitaemia=xs,site = -1, country= -1, study= -1, lambda_para=lambda_para),type='response', allow.new.levels = T)
plot(log10(xs), preds, type='l',
ylab = 'Pigment containing monocytes >0% (probability)',
xlab = 'Parasite density (per uL)', xlim = c(1.5,6.2),
xaxt='n', panel.first=grid(),lwd=3)
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
7 Pigment and severe anaemia
pigmt_data$SMA = as.numeric(pigmt_data$hb<5)
table(SMA=pigmt_data$SMA, study=pigmt_data$study)## study
## SMA AQ AQUAMAT Lyke et al SEAQUAMAT SMAC
## 0 457 2688 142 1168 20050
## 1 26 1158 24 100 5451par(las=1, mfrow=c(2,2))
plot(jitter(log10(pigmt_data$prop_pigneut[Af_ind]+lambda_pigment),amount = 0.05),xaxt='n',
jitter(pigmt_data$hb[Af_ind],amount = 0.02),
pch='.',ylim = c(2,13),
ylab='Haemoglobin (g/dL)', xlab = 'Pigment contaning PMNs (%)')
axis(1, at = log10(c(0,1,2,5,10,20)/100 +lambda_pigment),
labels = c(0,1,2,5,10,20))
mod = gam(hb ~ s(log10(prop_pigneut+lambda_pigment),k=5),
data = pigmt_data[Af_ind, ])
summary(mod)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hb ~ s(log10(prop_pigneut + lambda_pigment), k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.36439 0.01557 473 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(log10(prop_pigneut + lambda_pigment)) 2.898 3.307 242.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0266 Deviance explained = 2.67%
## GCV = 7.1549 Scale est. = 7.1539 n = 29510xs = (0: round(100*quantile(pigmt_data$prop_pigneut,
probs = .99, na.rm = T)))/100
lines(log10(xs+lambda_pigment),
predict(mod, newdata=data.frame(prop_pigneut=xs,
lambda_pigment=lambda_pigment)),
lwd=3,col='red')
plot(jitter(log10(pigmt_data$prop_pigmono[Af_ind]+lambda_pigment),amount = 0.05),xaxt='n',
jitter(pigmt_data$hb[Af_ind],amount = 0.02),
pch='.',ylim = c(2,13),
ylab='Haemoglobin (g/dL)', xlab = 'Pigment contaning monocytes (%)')
axis(1, at = log10(c(0,1,5,20,40)/100+lambda_pigment),
labels = c(0,1,5,20,40))
mod = gam(hb ~ s(log10(prop_pigmono+lambda_pigment),k=5)+
s(country_f,bs='re'),
data = pigmt_data[Af_ind, ])
summary(mod)##
## Family: gaussian
## Link function: identity
##
## Formula:
## hb ~ s(log10(prop_pigmono + lambda_pigment), k = 5) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.3482 0.2332 31.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(log10(prop_pigmono + lambda_pigment)) 3.669 3.913 504.1 <2e-16 ***
## s(country_f) 11.753 12.000 146.0 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.145 Deviance explained = 14.6%
## GCV = 6.3152 Scale est. = 6.3116 n = 28668xs = (0: round(100*quantile(pigmt_data$prop_pigmono,
probs = .99, na.rm = T)))/100
lines(log10(xs+lambda_pigment),
predict(mod, newdata=data.frame(prop_pigmono=xs,
lambda_pigment=lambda_pigment,
country_f= -1),
exclude = 'country_f'),lwd=3,col='red')## Warning in predict.gam(mod, newdata = data.frame(prop_pigmono = xs,
## lambda_pigment = lambda_pigment, : factor levels -1 not in original fit# SMA
mod = gam(SMA ~ s(log10(prop_pigneut+lambda_pigment),k=3)+
s(country_f, bs='re'),
family = 'binomial',
data = pigmt_data[Af_ind, ])
summary(mod)##
## Family: binomial
## Link function: logit
##
## Formula:
## SMA ~ s(log10(prop_pigneut + lambda_pigment), k = 3) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1394 0.1624 -7.017 2.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(prop_pigneut + lambda_pigment)) 1.961 1.998 43.57 <2e-16 ***
## s(country_f) 11.582 12.000 1124.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0441 Deviance explained = 4.38%
## UBRE = 0.019976 Scale est. = 1 n = 29510xs = (0: round(100*quantile(pigmt_data$prop_pigneut[Af_ind],
probs = .99, na.rm = T)))/100
preds = predict(mod, newdata=data.frame(prop_pigneut=xs,
lambda_pigment=lambda_pigment,
country_f= -1),
type='response',se.fit = T,exclude='country_f')## Warning in predict.gam(mod, newdata = data.frame(prop_pigneut = xs,
## lambda_pigment = lambda_pigment, : factor levels -1 not in original fitplot(log10(xs+lambda_pigment),
preds$fit,lwd=3,col='red',type='l',xaxt='n')
axis(1, at = log10(c(0,1,2,5,10,20)/100 +lambda_pigment),
labels = c(0,1,2,5,10,20))
mod = gam(SMA ~ s(log10(prop_pigmono+lambda_pigment),k=3)+
s(country_f, bs='re'),
family = 'binomial',
data = pigmt_data[Af_ind, ])
summary(mod)##
## Family: binomial
## Link function: logit
##
## Formula:
## SMA ~ s(log10(prop_pigmono + lambda_pigment), k = 3) + s(country_f,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2299 0.1749 -7.03 2.06e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(prop_pigmono + lambda_pigment)) 1.98 2 530.3 <2e-16 ***
## s(country_f) 11.62 12 841.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0586 Deviance explained = 6.11%
## UBRE = -0.0028473 Scale est. = 1 n = 28668xs = (0: round(100*quantile(pigmt_data$prop_pigmono[Af_ind],
probs = .99, na.rm = T)))/100
preds = predict(mod, newdata=data.frame(prop_pigmono=xs,lambda_pigment=lambda_pigment,
country_f= -1),
type='response',se.fit = T,exclude='country_f')## Warning in predict.gam(mod, newdata = data.frame(prop_pigmono = xs,
## lambda_pigment = lambda_pigment, : factor levels -1 not in original fitplot(log10(xs+lambda_pigment),
preds$fit,lwd=3,col='red',type='l',xaxt='n')
axis(1, at = log10(c(0,1,2,5,10,20)/100 +lambda_pigment),
labels = c(0,1,2,5,10,20))
8 HRP2 analysis - AQUAMAT
This can only be run with access to the HRP2 data
hrp2_aq = read.csv('~/Dropbox/Datasets/AQUAMAT/AQM_childHRP2.csv')
pigmt_data = pigmt_data %>% filter(study == 'AQUAMAT')
pigmt_data = merge(pigmt_data, hrp2_aq, by.x = 'id', by.y = 'studynumber')pigmt_data =
pigmt_data %>% filter(!is.na(hrp2_plasma),
!is.na(pig_neut) |
!is.na(pig_mono) |
!is.na(parasitaemia))
pigmt_data$SM = as.numeric(pigmt_data$hrp2_plasma >= 1000)
cor.test(pigmt_data$SM, as.numeric(pigmt_data$prop_pigneut>0.05))##
## Pearson's product-moment correlation
##
## data: x and y
## t = 8.6382, df = 2931, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1220675 0.1926587
## sample estimates:
## cor
## 0.1575644table(pigmt_data$SM)##
## 0 1
## 1049 1890sum(!is.na(pigmt_data$hrp2_plasma) & !is.na(pigmt_data$pig_neut))## [1] 2933sum(!is.na(pigmt_data$hrp2_plasma) & !is.na(pigmt_data$pig_mono))## [1] 2913mod=glmer(outcome ~ (prop_pigneut>0.05) * SM + (1 | site), family='binomial', data = pigmt_data)
summary(mod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: outcome ~ (prop_pigneut > 0.05) * SM + (1 | site)
## Data: pigmt_data
##
## AIC BIC logLik deviance df.resid
## 1807.6 1837.5 -898.8 1797.6 2928
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5097 -0.3519 -0.2957 -0.2206 5.3942
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.1525 0.3905
## Number of obs: 2933, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.78258 0.20103 -13.842 <2e-16 ***
## prop_pigneut > 0.05TRUE 0.68120 0.27249 2.500 0.0124 *
## SM 0.34810 0.17561 1.982 0.0475 *
## prop_pigneut > 0.05TRUE:SM -0.08223 0.30929 -0.266 0.7903
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pr_>0.05TRUE SM
## pr_>0.05TRUE -0.370
## SM -0.573 0.441
## p_>0.05TRUE: 0.328 -0.871 -0.559table(SM=pigmt_data$SM, PMNs_5percent=pigmt_data$PMNs5)## PMNs_5percent
## SM 0 1
## 0 872 172
## 1 1306 583table(SM=pigmt_data$SM, PMMs_5percent=pigmt_data$PMMs5)## PMMs_5percent
## SM 0 1
## 0 536 501
## 1 672 1204mod_para_hrp2 = gam(SM ~ s(log10(parasitaemia+50),k=5) + s(site_f, bs='re'), family='binomial', data = pigmt_data)
mod_PMN_hrp2 = glmer(SM ~ log10(prop_pigneut+lambda_pigment) + (1|site_f), family='binomial', data = pigmt_data)
mod_PMM_hrp2 = glmer(SM ~ log10(prop_pigmono+lambda_pigment) + (1|site_f), family='binomial', data = pigmt_data)
summary(mod_para_hrp2)##
## Family: binomial
## Link function: logit
##
## Formula:
## SM ~ s(log10(parasitaemia + 50), k = 5) + s(site_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4486 0.2459 1.825 0.0681 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(log10(parasitaemia + 50)) 2.805 3.32 52.86 <2e-16 ***
## s(site_f) 7.547 8.00 143.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0736 Deviance explained = 6.04%
## UBRE = 0.237 Scale est. = 1 n = 2798summary(mod_PMN_hrp2)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.8407528 0.22515897 8.175348 2.950141e-16
## log10(prop_pigneut + lambda_pigment) 0.8525308 0.08706886 9.791455 1.225199e-22summary(mod_PMM_hrp2)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1104826 0.19303018 5.752896 8.772726e-09
## log10(prop_pigmono + lambda_pigment) 0.5039785 0.06196779 8.132912 4.190991e-16Nsims=100
xs = (0:40)/100
preds_PMN =
bootMer(mod_PMN_hrp2, nsim=Nsims,
FUN=function(x) {
predict(x, newdata=data.frame(prop_pigneut=xs,
lambda_pigment=lambda_pigment),
re.form=NA)
})
preds_PMM =
bootMer(mod_PMM_hrp2, nsim=Nsims,
FUN=function(x){
predict(x, newdata=data.frame(prop_pigmono=xs,
lambda_pigment=lambda_pigment),
re.form=NA)
})ymax = quantile(pigmt_data$hrp2, probs = .995)
par(las = 1, mfrow=c(2, 3), mar = c(5,5,3,2),
family='serif',cex.lab=1.3, cex.axis=1.3)
plot(log10(pigmt_data$parasitaemia+50), log10(pigmt_data$hrp2+1),
col = adjustcolor(c('blue','red')[pigmt_data$outcome+1],
alpha.f = c(.2)),
pch=20, panel.first=grid(NA,NULL),ylim = log10(c(1, ymax)),
yaxt = 'n', xaxt='n', ylab='PfHRP2 (ng/mL)',
xlab = 'Parasite count (per uL)')
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
abline(h = log10(1000), v = 5, lty=2)
axis(2, at = 0:4, labels = c(1,10,
expression(10^2),
expression(10^3),
expression(10^4)))
mtext(text = 'a',side = 3,adj = 0, cex = 1.2, line=1.1)
xs = jitter(log10(pigmt_data$prop_pigneut+lambda_pigment*2),amount = .02)
plot(xs,log10(pigmt_data$hrp2+1), xlim = c(-2.05, log10(0.4+lambda_pigment*2)),
col = adjustcolor(c('blue','red')[pigmt_data$outcome+1],
alpha.f = c(.2)),
pch=20, panel.first=grid(NA,NULL),ylim = log10(c(1, ymax)),
yaxt = 'n', xaxt='n', ylab='PfHRP2 (ng/mL)',
xlab = 'Pigment containing PMNs (%)')
axis(1, at = log10(c(0,1,5,10,20,40)/100 + lambda_pigment*2), labels = c(0,1,5,10,20,40))
abline(v=log10(c(0,1,5,10,20,40)/100+lambda_pigment*2), lty='dotted', col='lightgray')
axis(2, at = 0:4, labels = c(1,10,
expression(10^2),
expression(10^3),
expression(10^4)))
legend('bottomright', fill=adjustcolor(c('blue','red'),.5), inset = 0.03,
legend = c('survived','died'))
abline(h = log10(1000), v = log10(0.05+lambda_pigment*2), lty=2)
mtext(text = 'b',side = 3,adj = 0, cex = 1.2, line=1.1)
xs = jitter(log10(pigmt_data$prop_pigmono+lambda_pigment*10),amount = .02)
plot(xs, log10(pigmt_data$hrp2+1), xlim = c(-1.32, log10(0.4+lambda_pigment*10)),
col = adjustcolor(c('blue','red')[pigmt_data$outcome+1],
alpha.f = c(.2)),
pch=20, panel.first=grid(NA,NULL),ylim = log10(c(1, ymax)),
yaxt = 'n', xaxt='n', ylab='PfHRP2 (ng/mL)',
xlab = 'Pigment containing monocytes (%)')
axis(1, at = log10(c(0,5,10,20,40)/100+lambda_pigment*10), labels = c(0,5,10,20,40))
abline(v=log10(c(0,5,10,20,40)/100+lambda_pigment*10), lty='dotted', col='lightgray')
axis(2, at = 0:4, labels = c(1,10,
expression(10^2),
expression(10^3),
expression(10^4)))
abline(h = log10(1000), v = log10(0.05+lambda_pigment*10), lty=2)
mtext(text = 'c',side = 3,adj = 0, cex = 1.2, line=1.1)
xs=seq(2,6,length.out=100)
my_preds = predict(mod_para_hrp2, newdata = data.frame(parasitaemia = 10^(xs), site_f = 'Kilifi'), se.fit=T)
plot(xs, 100*inv.logit(my_preds$fit),
lwd=3, xaxt='n', xlab = 'Parasite count (per uL))',
ylab = 'PfHRP2>1000 ng/mL (%)', panel.first = grid(),
ylim = c(0, 100),type='l')
axis(1, at = 2:6, labels = c(expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6)))
polygon(c(xs, rev(xs)),
100*inv.logit(c(my_preds$fit+1.96*my_preds$se.fit,
rev(my_preds$fit-1.96*my_preds$se.fit))),
col = adjustcolor('grey',.3), border = NA)
mtext(text = 'd',side = 3,adj = 0, cex = 1.2, line=1.1)
xs = (0:40)/100
plot(log10(xs+lambda_pigment*2), 100*inv.logit(apply(preds_PMN$t,2,median)),
lwd=3, xaxt='n', xlab = 'Pigment containing PMNs (%)',
ylab = 'PfHRP2>1000 ng/mL (%)', panel.first = grid(NA,NULL),
ylim = c(0, 100),type='l')
axis(1, at = log10(c(0,1,5,10,20,40)/100 +lambda_pigment*2),
labels = c(0,1,5,10,20,40))
polygon(log10(c(xs, rev(xs))+lambda_pigment*2),
100*inv.logit(c(apply(preds_PMN$t,2,quantile,probs=0.975),
rev(apply(preds_PMN$t,2,quantile,probs=0.025)))),
col = adjustcolor('grey',.3), border = NA)
abline(v=log10(c(0,1,5,10,20,40)/100+lambda_pigment*2), lty='dotted', col='lightgray')
mtext(text = 'e',side = 3,adj = 0, cex = 1.2, line=1.1)
plot(log10(xs+lambda_pigment*10),100*inv.logit(apply(preds_PMM$t,2,median)),
lwd=3, xaxt='n', xlab = 'Pigment containing monocytes (%)',
ylab = 'PfHRP2>1000 ng/mL (%)', panel.first = grid(NA,NULL),
ylim = c(0, 100),type='l')
axis(1, at = log10(c(0,5,10,20,40)/100 +lambda_pigment*10),
labels = c(0,5,10,20,40))
polygon(log10(c(xs, rev(xs))+lambda_pigment*10),
100*inv.logit(c(apply(preds_PMM$t,2,quantile,probs=0.975),
rev(apply(preds_PMM$t,2,quantile,probs=0.025)))),
col = adjustcolor('grey',.3), border = NA)
abline(v=log10(c(0,5,10,20,40)/100+lambda_pigment*10), lty='dotted', col='lightgray')
mtext(text = 'f',side = 3,adj = 0, cex = 1.2, line=1.1)